Vignette of analysis: Early-life demographic processes do not drive adult sex ratio biases and mating systems in sympatric coucal species

Supplementary Material A

Authors
Affiliations

Behavioural Genetics and Evolutionary Ecology, Max Planck Institute for Biological Intelligence

Department of Ornithology, Max Planck Institute for Biological Intelligence

Ignas Safari

Department Behavioural Neurobiology, Max Planck Institute for Biological Intelligence

Coucal Project

Department of Biology, College of Natural and Mathematical Sciences, University of Dodoma

Poyo Makomba

Coucal Project

Anne Hertel

Department Biology II, Ludwig Maximilians University Munich

Department Behavioural Neurobiology, Max Planck Institute for Biological Intelligence

Coucal Project

Department Biology II, Ludwig Maximilians University Munich

Published

November 22, 2023

In this document we provide all the necessary code for reproducing the analyses presented in our manuscript. To access the dataset and Rmarkdown file, please download this GitHub repository. Simply follow the link and click on Download ZIP on the right-hand side of the page. An explanation of the files in the repository can be found in the Readme file. Please don’t hesitate to contact Luke at luke.eberhart[at]bi.mpg.de if you have any questions.

Prerequisites

GitHub repository

  • For running the complete code you need a files subfolder containing the raw data downloaded from Data_files folder provided in the GitHub repository.

R packages

  • The following packages are needed for analysis and can be easily installed from CRAN using the following code:
# a vector of all the packages needed in the project
packages_required_in_project <- c("RMark",
                                  "tidyverse",
                                  "readxl",
                                  "BaSTA",
                                  "pbapply",
                                  "RColorBrewer",
                                  "grid",
                                  "Rmisc",
                                  "gss",
                                  "arm",
                                  "partR2",
                                  "parameters",
                                  "standardize",
                                  "colorBlindness",
                                  "ggthemes",
                                  "patchwork",
                                  "gt",
                                  "rptR",
                                  "tidybayes",
                                  "broom.mixed",
                                  "effects",
                                  "patchwork",
                                  "devtools",
                                  "unmarked",
                                  "R2ucare",
                                  "marked",
                                  "merTools",
                                  "bootpredictlme4",
                                  "extrafont",
                                  "survminer",
                                  "bdsmatrix",
                                  "coxme",
                                  "epitools",
                                  "survival",
                                  "magrittr",
                                  "ggpubr",
                                  "reshape2",
                                  "sp",
                                  "adehabitatLT",
                                  "reshape",
                                  "coefplot",
                                  "mapview",
                                  "lubridate",
                                  "effects",
                                  "merDeriv",
                                  "remotes")
                                  
# of the required packages, check if some need to be installed
new.packages <- 
  packages_required_in_project[!(packages_required_in_project %in% 
                                   installed.packages()[,"Package"])]

# install all packages that are not locally available
if(length(new.packages)) install.packages(new.packages)

# install the bootpredictlme4 package directly from GitHub
remotes::install_github("RemkoDuursma/bootpredictlme4")

# load all the packages into the current R session
lapply(packages_required_in_project, require, character.only = TRUE)

Plotting themes

  • The following plotting themes, colors, and typefaces are used throughout the project:
# Find fonts from computer that you want. Use regular expressions to do this
# For example, load all fonts that are 'candara' or 'Candara'
extrafont::font_import(pattern = "[V/v]erdana", prompt = FALSE)

# check which fonts were loaded
extrafont::fonts()
extrafont::fonttable()
extrafont::loadfonts() # load these into R

# define the plotting theme to be used in subsequent ggplots
luke_theme <- 
  theme_bw() +
  theme(
    text = element_text(family = "Verdana"),
    legend.title = element_text(size = 16),
    legend.text = element_text(size = 12),
    axis.title.x = element_text(size = 12),
    axis.text.x  = element_text(size = 10), 
    axis.title.y = element_text(size = 12),
    axis.text.y = element_text(size = 10),
    strip.text = element_text(size = 12),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.ticks = element_line(size = 0.5, colour = "grey40"),
    axis.ticks.length = unit(0.2, "cm"),
    panel.border = element_rect(linetype = "solid", colour = "grey"),
    legend.position = c(0.1, 0.9)
  )

# set plotting color palettes
sex_pal2 <- 
  c(pull(ggthemes_data$wsj$palettes$colors6[3,2]), 
    pull(ggthemes_data$wsj$palettes$colors6[2,2]))
sex_pal3 <- 
  c(pull(ggthemes_data$wsj$palettes$colors6[3,2]), 
    pull(ggthemes_data$wsj$palettes$colors6[3,2]),
    pull(ggthemes_data$wsj$palettes$colors6[2,2]),
    pull(ggthemes_data$wsj$palettes$colors6[2,2]))
plot_palette_sex <- RColorBrewer::brewer.pal(8, "Dark2")[c(2,1)]
plot_palette_species <- RColorBrewer::brewer.pal(8, "Dark2")[c(3,7)]

# specify the facet labels for each species and sex
species_names <- c(
  'BC' = "Black coucal",
  'WBC' = "White-browed coucal")
sex_names <- c(
  'female' = "Females",
  'male' = "Males")

# color of mean estimate point in forest plots
col_all <- "#2E3440"

Analytical Functions

The following custom functions are used to bootstrap the various datasets and run the stochastic demographic model. To reproduce the results presented in our manuscript, these functions need to be stored in the environment of your R Session.

adult_CJS_bootstrap_data()

This function randomly samples the adult mark-recapture dataset to produce a random one the same size as the original, with replacement

arguments:

  • adult: the adult mark-recapture dataset
  • num_boot: a unique number for the pbapply simulation (don’t specify)
  • species: species name (either “BC” or “WBC”)
  • iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
adult_CJS_bootstrap_data <- 
  function(adult, num_boot, species, iter_add) {
  
  # sample a new adult mark-recapture dataset the same size as the original, 
  # with replacement
  adult_boot <-
    adult %>%
    sample_n(nrow(.), replace = TRUE)
  
  # make a list of the new dataset, iteration, species
  out <- list(adult_boot = adult_boot,
              iter = num_boot + ((iter_add - 1) * niter),
              species = species)
}

bootstrap_adult_survival_CJS()

runs the CJS survival analyses using the the bootstrapped sample created from adult_CJS_bootstrap_data().

arguments:

  • coucal_boot_list: the bootstrapped adult mark-recapture dataset produced from the adult_CJS_bootstrap_data() function above
  • num_boot: a unique number for the pbapply simulation (don’t specify)
  • first_year: first year of the study (i.e., the year of the first event in the encounter history)
  • bootstrap_name: a unique string for naming of files in the output directory
  • species: species name (either “BC” or “WBC”)
  • iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
  • prefix_number: a unique prefix to prevent overwriting output data on disk (e.g., “boot_one_”)
  • out_directory: location in the project directory where the output files should be stored (e.g., “tmp”)
bootstrap_adult_survival_CJS <- 
  function(coucal_boot_list, 
           num_boot, 
           first_year,
           bootstrap_name,
           species,
           iter_add,
           prefix_number,
           out_directory = "/output/bootstraps/raw/") {
    
    # specify the bootstrapped data samples (from the previous function)
    adult_ch <- coucal_boot_list[["adult_boot"]]
    
    # clean up capture histories
    adult_ch <-    
      adult_ch %>% 
      ungroup() %>% 
      as.data.frame() %>% 
      mutate(sex = as.factor(sex))
    
    # process the adult data as a "CJS"  analysis
    coucal_adult.proc <- RMark::process.data(adult_ch, 
                                             model = "CJS",
                                             groups = c("sex"),
                                             begin.time = first_year)
    
    # Create the design matrix from the processed mark-recapture datasets
    coucal_adult.ddl <- RMark::make.design.data(coucal_adult.proc)
    
    adult_survival_analysis_run = 
      function(proc_data, design_data){
        # sex- and stage-specific survival:
        Phi.sex = list(formula = ~ sex) 
        
        # Models exploring variation in encounter probability
        # constant:
        p.dot = list(formula =  ~ 1)
        
        # sex-dependent:
        p.sex = list(formula =  ~ sex)
        
        # factorial variation across year:
        p.Time = list(formula =  ~ time)
        
        # additive effects of sex and factorial year:
        p.sex_Time = list(formula =  ~ sex + time)
        
        # create a list of candidate models for all the a models above that begin with 
        # either "Phi." or "p."
        cml <-  RMark::create.model.list("CJS")
        
        # specify the data, design matrix, delete unneeded output files, and 
        # run the models in Program MARK
        model.list <-  RMark::mark.wrapper(cml, data = proc_data, 
                                           ddl = design_data, delete = TRUE, 
                                           wrap = FALSE, threads = 1, brief = TRUE,
                                           silent = TRUE, output = FALSE, prefix = prefix_number)
        
        # output the model list and sotre the results
        return(model.list)
      }
    
    # Run the models on the bootstrapped data
    adult_survival_analysis_out <-
      adult_survival_analysis_run(proc_data = coucal_adult.proc,
                                  design_data = coucal_adult.ddl)
    
    extract_top_model_output <- 
      function(rmark_output, top_model = TRUE, mod_num){
        # Find the model number for the first ranked model of the AIC table
        if(top_model == TRUE){
          mod_num <- 
            as.numeric(rownames(rmark_output$model.table[1,]))
        }
        
        else{
          mod_num <- mod_num
        }
        
        # extract and wrangle reals from model output 
        reals <- 
            rmark_output[[mod_num]]$results$real %>% 
            bind_cols(data.frame(str_split_fixed(rownames(.), " ", n = 5)), .) %>% 
            filter(X1 == "Phi") %>% 
            mutate(sex = as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) == "F", "Female", "Male"))) %>% 
            dplyr::select(sex, estimate) %>% 
            mutate(iter = num_boot + ((iter_add - 1) * niter),
                   species = species)
        
        # extract and wrangle the beta estimates from the model output
        betas <- 
          rmark_output[[mod_num]]$results$beta %>% 
          mutate(statistic = rownames(.),
                 iter = num_boot + ((iter_add - 1) * niter),
                 species = species) %>% 
          mutate(parameter = ifelse(grepl(x = statistic, pattern = "S:"), "S",
                                    ifelse(grepl(x = statistic, pattern = "p:"), "p",
                                           ifelse(grepl(x = statistic, pattern = "r:"), "r",
                                                  ifelse(grepl(x = statistic, pattern = "F:"), "F", 
                                                         ifelse(grepl(x = statistic, pattern = "Phi:"), "Phi", "XXX"))))),
                 variable = ifelse(grepl(x = statistic, pattern = "Intercept"), "Intercept",
                                   ifelse(grepl(x = statistic, pattern = "sexM:Cubic"), "sexM:Cubic",
                                          ifelse(grepl(x = statistic, pattern = "sexM:Quadratic"), "sexM:Quadratic",
                                                 ifelse(grepl(x = statistic, pattern = "sexM:Time"), "sexM:Time",
                                                        ifelse(grepl(x = statistic, pattern = "sexM"), "sexM",
                                                               ifelse(grepl(x = statistic, pattern = "Time"), "Time",
                                                                      ifelse(grepl(x = statistic, pattern = "Cubic"), "Cubic", 
                                                                             ifelse(grepl(x = statistic, pattern = "Quadratic"), "Quadratic","XXX"))))))))) %>% 
          dplyr::select(-statistic)
        
        # consolidate the AIC model selection results
        AIC_table <- 
          rmark_output$model.table %>% 
          mutate(iter = num_boot + ((iter_add - 1) * niter),
                 species = species,
                 model_no_orig = as.numeric(rownames(.))) %>% 
          mutate(model_no_rank = as.numeric(rownames(.)))
        
        # consolidate the all output into a list
        survival_model_output_list <- 
          list(reals = reals,
               betas = betas,
               AIC_table = AIC_table)
        
        survival_model_output_list
      }
    
    # extract and format survival rates from model output
    adult_survival_model_output_list <- 
      extract_top_model_output(rmark_output = adult_survival_analysis_out)
    
    # make a list of all the results from this iteration
    bootstrap_results_list <- 
      list(adult_survival_model_output = adult_survival_model_output_list,
           bootstrapped_data = coucal_boot_list)
    
    # extract file directory
    file_directory <- paste0(getwd(), "/", out_directory, "/", bootstrap_name, "_", num_boot, ".Rds")
    
    # save the results of the iteration
    save(bootstrap_results_list, file = file_directory)
    
    bootstrap_results_list
  }

run_bootstrap_adult_survival_CJS

run the bootstrap and the CJS analysis (passed through the pbsapply function), see arguments defined abov

arguments:

  • adult: the adult mark-recapture dataset
  • num_boot: a unique number for the pbapply simulation (don’t specify)
  • species: species name (either “BC” or “WBC”)
  • iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
  • first_year: first year of the study (i.e., the year of the first event in the encounter history)
  • bootstrap_name: a unique string for naming of files in the output directory
  • prefix_number: a unique prefix to prevent overwriting output data on disk (e.g., “boot_one_”)
  • out_directory: location in the project directory where the output files should be stored (e.g., “tmp”)
run_bootstrap_adult_survival_CJS <- function(adult, num_boot, species, iter_add,
                                     first_year, bootstrap_name, prefix_number,
                                     out_directory)
{
  # run the sampling function and specify the datasets
  bootstrap_data_list <- 
    adult_CJS_bootstrap_data(adult = adult, 
                             num_boot = num_boot,
                             species = species, 
                             iter_add = iter_add)
  
  # run the survival analysis and ASR deduction on the sampled data
  result <- bootstrap_adult_survival_CJS(coucal_boot_list = bootstrap_data_list, 
                                         num_boot = num_boot, 
                                         first_year = first_year,
                                         bootstrap_name = bootstrap_name,
                                         species = species,
                                         iter_add = iter_add,
                                         prefix_number = prefix_number,
                                         out_directory = out_directory)
  result
}

surv_boot_out_wrangle()

loads all the output rds files from the parallel adult CJS bootstrapping procedure and consolidates them into a single object for further analysis

arguments:

  • species: species name (either “BC” or “WBC”)
  • niter: number of iterations specified in the stochastic simulation (used to predict how many files to merge from the output directory)
  • out_directory: location in the project directory where the output files from the simulation were stored (e.g., “tmp”)
# boot_out_wrangle() 
surv_boot_out_wrangle <- function(species, niter, output_directory = "output/bootstraps/"){
  
  # specify directory string
  boot_out_path <- paste0(output_directory, 
                          species,"_adult_survival_CJS_bootstrap_result.rds")
  
  # load the bootstrap output file
  survival_bootstrap_result <- readRDS(boot_out_path)
  
  # bind all the adult rates output together
  adult_reals_survival_rates_boot <-
    rbind(do.call(rbind, lapply(seq(from = 1, to = niter, by = 1),
                                function(x) survival_bootstrap_result["adult_survival_model_output", x]$adult_survival_model_output$reals))) %>%
    arrange(iter, sex)
  
  # bind all the adult beta estimates output together
  adult_betas_survival_rates_boot <-
    rbind(do.call(rbind, lapply(seq(from = 1, to = niter, by = 1), 
                                function(x) survival_bootstrap_result["adult_survival_model_output", x]$adult_survival_model_output$betas)))
  
  # bind all the adult AIC table output together
  adult_AIC_tables_boot <-
    rbind(do.call(rbind, lapply(seq(from = 1, to = niter, by = 1), 
                                function(x) survival_bootstrap_result["adult_survival_model_output", x]$adult_survival_model_output$AIC_table)))
  
  # tidy all the bound data together into a mega list of results
  boot_out_tidy <- 
    list(adult_reals_survival_rates_boot = adult_reals_survival_rates_boot,
         adult_betas_survival_rates_boot = adult_betas_survival_rates_boot,
         adult_AIC_tables_boot = adult_AIC_tables_boot)
  
  boot_out_tidy
}

approx_sd()

derive the approximate sd given the upper and lower 95% confidence intervals (see:Cochrane Chapter 7.7.3.2)

arguments:

  • x1: lower bound of the confidence interval
  • x2: upper bound of the confidence interval
# function to approximate the sd given the 95% CIs
approx_sd <- function(x1, x2){
  (x2-x1) / (qnorm(0.975) - qnorm(0.025))
}

matrix_ASR()

calculates the ASR of the population based on the two-sex two-stage projection matrix built by the coucal_matrix() function (below).

arguments:

  • M: a two-sex x by x projection matrix produced by the coucal_matrix() function above
  • n: x-lengthed vector representing starting stage distribution (the default is a vector with 10 individuals in each stage)
  • h: harem size index (default is 1, monogamy; polyandry < 1, polygyny > 1)
  • k: clutch size (default is 4 eggs, mean for both coucal species)
  • iterations: number of time steps to project the demogrphic model (n.b., after 100 steps our model reaches the stable stage distribution)
  • HSR: hatching sex ratio (expressed as the proportion of male chicks)
  • num_boot: a unique number for the pbapply simulation (don’t specify)
  • species: species name (either “BC” or “WBC”)
  • iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk (if you run this on one session, then don’t specify)
matrix_ASR <-
  function(M, 
           n = rep(10, nrow(M)), 
           h = 1, 
           k = 4, 
           iterations = 1000, 
           HSR = 0.5, 
           num_boot, 
           species, 
           iter_add){
    
    # Number of stages in matrix
    x <- length(n)
    
    # Number of time steps to simulate
    t <- iterations
    
    # an empty t by x matrix to store the stage distributions
    stage <- matrix(data = numeric(x * t), nrow = x, ncol = t)
    
    # an empty t by x matrix to store the stage-specific frequencies
    pop_dist <- matrix(data = numeric(x * t), nrow = x, ncol = t)
    rownames(pop_dist) <- rownames(M)
    colnames(pop_dist) <- 0:(t - 1)
    
    # an empty t vector to store the population sizes
    pop <- numeric(t)
    
    # for loop that goes through each of t time steps
    for (i in 1:t) {
      
      # stage distribution at time t
      stage[,i] <- n
      
      # population size at time t
      pop[i] <- sum(n)
      
      # number of male adults at time t = (number alive at t-1) * (survival rate)
      M2 <- stage[4, i] * M["M_Adult", "M_Adult"]
      
      # number of female adults at time t
      F2 <- stage[2, i] * M["F_Adult", "F_Adult"]
      
      # Female freq-dep fecundity of Female chicks
      M["F_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * (1 - HSR))
      
      # Female freq-dep fecundity of Male chicks
      M["M_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * HSR)
      
      # Male freq-dep fecundity of Female chicks
      M["F_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * (1 - HSR))
      
      # Male freq-dep fecundity of Male chicks
      M["M_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * HSR)
      
      # define the new n (i.e., new stage distribution at time t)
      n <- M %*% n
      
      # define rownames of stage matrix
      rownames(stage) <- rownames(M)
      
      # define colnames of stage matrix
      colnames(stage) <- 0:(t - 1)
      
      # calculate the proportional stable stage distribution
      stage <- apply(stage, 2, function(x) x/sum(x))
      
      # define stable stage as the last stage
      stable.stage <- stage[, t]
      
      pop_dist[, i] <- n
      
    }
    # calc ASR as the proportion of the adult stable stage class that is male
    ASR <- stable.stage[4]/(stable.stage[2] + stable.stage[4])
    SSR <- stable.stage[3]/(stable.stage[1] + stable.stage[3])
    
    # make a list of results
    pop.proj <- list(ASR = ASR,
                     SSR = SSR,
                     lambda = pop[t]/pop[t - 1],
                     pop_size = pop,
                     stage_size = pop_dist,
                     stable.stage = stable.stage,
                     stage.vectors = stage,
                     SSD_M2 = stable.stage[4],
                     SSD_F2 = stable.stage[2],
                     iter = num_boot + ((iter_add - 1) * niter), 
                     species = species)
    
    # print the list as output to the function
    pop.proj
  }

coucal_matrix()

builds the two-sex Lefkovitch matrix using the vital rates specified in the demographic_rates object. This function also give the option to specify a one-sex matrix for analytical comparisons between 2- and 1-sex matrix models

arguments:

  • demographic_rates_: a list of all the vital rates of an iteration of the stochastic simulation run within the run_bootstrap_juv_hazd_ad_surv_ASR() function below
  • two_sex: whether to make a one-sex or two-sex matrix of the rates provided (default is ‘two_sex = TRUE’)
coucal_matrix <- 
  function(demographic_rates_, two_sex = TRUE){
    if(two_sex){
      
      # Define coucal life-stages of the coucal matrix model
      stages <- c("F_1st_year",  "F_Adult",  "M_1st_year",  "M_Adult")
      
      # Build the 4x4 matrix
      result <- 
        matrix(c(
          
          # top row of matrix
          0, NA, 0, NA, 
          
          # second row of matrix
          (demographic_rates_$Egg_survival * 
             demographic_rates_$F_Nestling_survival * 
             demographic_rates_$F_Groundling_survival *
             demographic_rates_$F_Fledgling_survival),
          demographic_rates_$F_Adult_survival,
          0, 0,
          
          # third row of matrix
          0, NA, 0, NA, 
          
          # fourth row of matrix
          0, 0, 
          (demographic_rates_$Egg_survival * 
             demographic_rates_$M_Nestling_survival * 
             demographic_rates_$M_Groundling_survival *
             demographic_rates_$M_Fledgling_survival),
          demographic_rates_$M_Adult_survival),
          nrow = length(stages), byrow = TRUE,
          dimnames = list(stages, stages))
    }
    else{
      
      # Define coucal life-stages of the Ceuta snowy coucal matrix model
      stages <- c("1st_year",  "Adult")
      
      # Build the 4x4 matrix
      result <- 
        matrix(c(
          
          # top row of matrix
          0, NA,
          
          # second row of matrix
          (demographic_rates_$Egg_survival * 
             demographic_rates_$Nestling_survival * 
             demographic_rates_$Groundling_survival *
             demographic_rates_$Fledgling_survival),
          demographic_rates_$Adult_survival),
          nrow = length(stages), byrow = TRUE,
          dimnames = list(stages, stages))
    }
    result
  }

bootstrap_hazard_data()

This function randomly samples one sibling per nest from the offspring survival dataset

arguments:

  • offspring: the radio-tracked offspring dataset
  • num_boot: a unique number for the pbapply simulation (don’t specify)
  • species: species name (either “BC” or “WBC”)
  • iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
  • alpha_value: Parameter defining cross-validation score for smoothing parameter selection in the gss::sshzd function.
  • max_time: maximum number of days since hatch to model survival (‘70’ in the case of our analysis)
  • niter: number of iterations specified in the stochastic simulation
# bootstrap_hazard_data() 
bootstrap_hazard_data <- 
  function(offspring, num_boot, species, iter_add, alpha_value, max_time, niter) {
    
    # set attempt to 0 at start of each loop
    attempt <- 0
    
    max_time <- max_time
    
    hzd_ss_try_M_lambda <- 1.1
    hzd_ss_try_F_lambda <- 1.1

    # store simulated estimates only if lambda is less than 1.1
    while( (hzd_ss_try_M_lambda > -1.1 | hzd_ss_try_F_lambda > 1.1) && attempt <= 100 ) {

      # next attempt
      attempt <- attempt + 1

      # sample a new offspring dataset containing only one nest member per draw
      offspring_boot <-
        offspring %>%
        group_by(nest_ID) %>%
        sample_n(1)

      try(
        hzd_ss_try_M <- sshzd(Surv(exit, event, entry) ~ exit,
                              data = filter(offspring_boot, sex == "M"),
                              alpha = alpha_value),
        silent = TRUE
      )

      try(
        hzd_ss_try_F <- sshzd(Surv(exit, event, entry) ~ exit,
                              data = filter(offspring_boot, sex == "F"),
                              alpha = alpha_value),
        silent = TRUE
      )

      # store calculated peak (i.e., the apex of the polynomial curve)
      try(
        hzd_ss_try_M_lambda <- hzd_ss_try_M$lambda
      )
      # store calculated peak (i.e., the apex of the polynomial curve)
      try(
        hzd_ss_try_F_lambda <- hzd_ss_try_F$lambda
      )
    }
    
    # store simulated estimates only if the hzdcurve can be estimated
    while( (exists("hzd_curve_try_M") == FALSE | exists("hzd_curve_try_F") == FALSE) && attempt <= max_time) {
      
      time_vector <- seq(0, max_time - attempt, 1)
      
      # next attempt
      attempt <- attempt + 1
      
      try(
        hzd_ss_try_M <- sshzd(Surv(exit, event, entry) ~ exit, 
                              data = filter(offspring_boot, sex == "M"), 
                              alpha = alpha_value),
        silent = TRUE
      )
      
      try(
        hzd_ss_try_F <- sshzd(Surv(exit, event, entry) ~ exit, 
                              data = filter(offspring_boot, sex == "F"), 
                              alpha = alpha_value),
        silent = TRUE
      )
      
      # simulate an estimate
      try(
        hzd_curve_try_M <- 
          hzdcurve.sshzd(object = hzd_ss_try_M, time = time_vector, se = TRUE),
        silent = TRUE
      )
      try(
        hzd_curve_try_F <- 
          hzdcurve.sshzd(object = hzd_ss_try_F, time = time_vector, se = TRUE),
        silent = TRUE
      )
      
    }
    
    # make a list including elements that will be used in the next function
    out <- list(offspring_boot = offspring_boot, 
                iter = num_boot + ((iter_add - 1) * niter),
                species = species,
                time_vector = time_vector,
                lambdaM = hzd_ss_try_M_lambda,
                lambdaF = hzd_ss_try_F_lambda)
  }

run_bootstrap_juv_hazd_ad_surv_ASR()

initiates the bootstrap_hazard_data() and runs the juvenile survival analysis to derive sex- and stage-specific survival rates. Then constructs the matrix model with the rates and projects this into the future to acquire the stable stage distribution and derive the ASR.

arguments:

  • num_boot: a unique number for the pbapply simulation (don’t specify)
  • offspring: the radio-tracked offspring dataset
  • k_dist: a two-number vector containing the mean and sd of clutch size (output from the analyses below)
  • HSR_dist: a two-number vector containing the mean and sd of hatching sex ratio (output from the analyses below)
  • h_dist: a two-number vector containing the mean and sd of the harem size index (output from the analyses below)
  • egg_surv_dist: a two-number vector containing the mean and sd of egg survival (output from the analyses below)
  • flight_age_distF: a two-number vector containing the mean and sd of female-specific flight age (output from the analyses below)
  • flight_age_distM: a two-number vector containing the mean and sd of male-specific flight age (output from the analyses below)
  • fledge_age_distF: a two-number vector containing the mean and sd of female-specific fledge age (output from the analyses below)
  • fledge_age_distM: a two-number vector containing the mean and sd of male-specific fledge age (output from the analyses below)
  • bootstrap_name: a unique string for naming of files in the output directory
  • adult_surival_boot_out: the bootstrapped adult encounter history data
  • species: species name (either “BC” or “WBC”)
  • iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
  • prefix_number: a unique prefix to prevent overwriting output data on disk (e.g., “boot_one_”)
  • alpha_value: Parameter defining cross-validation score for smoothing parameter selection in the gss::sshzd function (default is 1.4).
  • max_time: maximum number of days since hatch to model survival (‘70’ in the case of our analysis)
  • niter: number of iterations specified in the stochastic simulation
  • out.directory: location in the project directory where the output files should be stored (e.g., “tmp”)
run_bootstrap_juv_hazd_ad_surv_ASR <- function(num_boot, 
                                               offspring,
                                               k_dist, 
                                               HSR_dist, 
                                               h_dist, 
                                               egg_surv_dist, 
                                               flight_age_distF,
                                               flight_age_distM,
                                               fledge_age_distF, 
                                               fledge_age_distM,
                                               bootstrap_name, 
                                               adult_surival_boot_out,
                                               species, 
                                               iter_add, 
                                               prefix_number,
                                               alpha_value = 1.4,
                                               max_time = 70,
                                               niter,
                                               output.directory)
{
  # run the sampling function and specify the datasets
  coucal_boot_list <- 
    bootstrap_hazard_data(offspring = offspring, 
                          num_boot = num_boot,
                          species = species, 
                          iter_add = iter_add,
                          alpha_value = alpha_value,
                          max_time = max_time,
                          niter = niter)
  
  # run the survival analysis and ASR deduction on the sampled data
  # specify the bootstrapped data samples (from the previous function)
    offspring_data <- coucal_boot_list[["offspring_boot"]]
    
    # clean up capture histories
    offspring_data <-    
      offspring_data %>% 
      ungroup() %>% 
      as.data.frame()
    
    # fit smoothed spline of hazard function for either sex
    M_haz_ss <- sshzd(Surv(exit, event, entry) ~ exit, 
                      data = filter(offspring_data, sex == "M"), 
                      alpha = alpha_value)

    F_haz_ss <- sshzd(Surv(exit, event, entry) ~ exit, 
                      data = filter(offspring_data, sex == "F"), 
                      alpha = alpha_value)

    haz_ss_lambda <- list(Male_haz_ss_lambda = M_haz_ss$lambda,
                          Female_haz_ss_lambda = F_haz_ss$lambda)
    
    # extract fitted estimates from the spline function
    M_haz_ss_curve <- 
      hzdcurve.sshzd(object = M_haz_ss, time = coucal_boot_list[["time_vector"]], se = TRUE)
    
    F_haz_ss_curve <- 
      hzdcurve.sshzd(object = F_haz_ss, time = coucal_boot_list[["time_vector"]], se = TRUE)
    
    haz_ss_curve <- 
      expand.grid(species = species, 
                  age = coucal_boot_list[["time_vector"]],
                  sex = c("Male", "Female")) %>% 
      mutate(fit = c(M_haz_ss_curve$fit, F_haz_ss_curve$fit),
             se = c(M_haz_ss_curve$se, F_haz_ss_curve$se)) %>% 
      mutate(estimate = 1 - fit,
             upper = 1 - fit * exp(1.96 * se),
             lower = 1 - fit / exp(1.96 * se),
             iter = num_boot)
    
    # transform the daily nestling survival (DCS) to apparent fledgling success
    # by calculating the product of all DCS estimates:
    
    fledge_ageM <- rnorm(1, fledge_age_distM[1], fledge_age_distM[2])
    fledge_ageF <- rnorm(1, fledge_age_distF[1], fledge_age_distF[2])
    
    flight_ageM <- rnorm(1, flight_age_distM[1], flight_age_distM[2])
    flight_ageF <- rnorm(1, flight_age_distF[1], flight_age_distF[2])
    
    coucal_nestling_survivalF <-
      haz_ss_curve %>% 
      filter(age <= fledge_ageF & sex == "Female") %>% 
      group_by(sex) %>%
      dplyr::summarise(value = prod(estimate), .groups = 'drop') %>%
      mutate(stage = "nestling",
             rate = "survival")
    
    coucal_nestling_survivalM <-
      haz_ss_curve %>% 
      filter(age <= fledge_ageM & sex == "Male") %>% 
      group_by(sex) %>%
      dplyr::summarise(value = prod(estimate), .groups = 'drop') %>%
      mutate(stage = "nestling",
             rate = "survival")
    
    coucal_nestling_survival <- 
      bind_rows(coucal_nestling_survivalF, coucal_nestling_survivalM)
    
    coucal_groundling_survivalF <-
      haz_ss_curve %>% 
      filter(age < (flight_ageF + fledge_ageF) & age > fledge_ageF & sex == "Female") %>% 
      group_by(sex) %>% 
      dplyr::summarise(value = prod(estimate), .groups = 'drop') %>% 
      mutate(stage = "groundling",
             rate = "survival")
    
    coucal_groundling_survivalM <-
      haz_ss_curve %>% 
      filter(age < (flight_ageM + fledge_ageM) & age > fledge_ageM & sex == "Male") %>% 
      group_by(sex) %>% 
      dplyr::summarise(value = prod(estimate), .groups = 'drop') %>% 
      mutate(stage = "groundling",
             rate = "survival")
    
    coucal_groundling_survival <- 
      bind_rows(coucal_groundling_survivalF, coucal_groundling_survivalM)
    
    coucal_fledgling_survivalF <-
      haz_ss_curve %>% 
      filter(age >= (flight_ageF + fledge_ageF) & sex == "Female") %>% 
      group_by(sex) %>% 
      dplyr::summarise(value = prod(estimate), .groups = 'drop') %>% 
      mutate(stage = "fledgling",
             rate = "survival")
    
    coucal_fledgling_survivalM <-
      haz_ss_curve %>% 
      filter(age >= (flight_ageM + fledge_ageM) & sex == "Male") %>% 
      group_by(sex) %>% 
      dplyr::summarise(value = prod(estimate), .groups = 'drop') %>% 
      mutate(stage = "fledgling",
             rate = "survival")
    
    coucal_fledgling_survival <- 
      bind_rows(coucal_fledgling_survivalF, coucal_fledgling_survivalM)
    
    coucal_adult_survival <-
      filter(adult_surival_boot_out, iter == sample(1:niter, 1)) %>%
      dplyr::select(-species, -iter) %>%
      mutate(stage = "adult",
             rate = "survival")
    
    egg_survival <- rnorm(1, egg_surv_dist[1], egg_surv_dist[2])
    
    coucal_egg_survival <- 
      data.frame(sex = NA,
                 value = egg_survival,
                 stage = c("egg"),
                 rate = c("survival"))
    
    coucal_mating_system <- 
      data.frame(sex = NA,
                 value = 1/rnorm(1, h_dist[1], h_dist[2]),
                 stage = c("h"),
                 rate = c("fecundity"))
    
    coucal_clutch_size <- 
      data.frame(sex = NA,
                 value = rnorm(1, k_dist[1], k_dist[2]),
                 stage = c("k"),
                 rate = c("fecundity"))
    
    coucal_HSR <- 
      data.frame(sex = NA,
                 value = rnorm(1, HSR_dist[1], HSR_dist[2]),
                 stage = c("HSR"),
                 rate = c("fecundity"))
    
    coucal_flight_age <- 
      data.frame(sex = c("Female", "Male"),
                 value = c(flight_ageF, flight_ageM),
                 stage = c("flight_age"),
                 rate = c("development"))
    
    coucal_fledge_age <- 
      data.frame(sex = c("Female", "Male"),
                 value = c(fledge_ageF, fledge_ageM),
                 stage = c("fledge_age"),
                 rate = c("development"))
    
    # Bind the juvenile and adult dataframe with the nestlings
    coucal_vital_rates <- 
      bind_rows(coucal_egg_survival,
                coucal_nestling_survival,
                coucal_groundling_survival,
                coucal_fledgling_survival,
                coucal_adult_survival,
                coucal_mating_system,
                coucal_clutch_size,
                coucal_HSR,
                coucal_flight_age,
                coucal_fledge_age) %>% 
      as.data.frame() %>% 
      mutate(iter = num_boot, #+ ((iter_add - 1) * niter),
             species = species) %>% 
      arrange(sex, rate, stage)
    
    # Create a list of demographic rates from the survival analyses above
    demographic_rates <- list(Egg_survival = pull(filter(coucal_vital_rates, stage == "egg"), value),
                              F_Nestling_survival = pull(filter(coucal_vital_rates, 
                                                                stage == "nestling" & rate == "survival" & sex == "Female"), 
                                                         value),
                              F_Groundling_survival = pull(filter(coucal_vital_rates, 
                                                                  stage == "groundling" & rate == "survival" & sex == "Female"), 
                                                           value),
                              F_Fledgling_survival = pull(filter(coucal_vital_rates, 
                                                                 stage == "fledgling" & rate == "survival" & sex == "Female"), 
                                                          value),
                              F_Adult_survival = pull(filter(coucal_vital_rates, 
                                                             stage == "adult" & rate == "survival" & sex == "Female"), 
                                                      value),
                              M_Nestling_survival = pull(filter(coucal_vital_rates, 
                                                                stage == "nestling" & rate == "survival" & sex == "Male"), 
                                                         value),
                              M_Groundling_survival = pull(filter(coucal_vital_rates, 
                                                                  stage == "groundling" & rate == "survival" & sex == "Male"), 
                                                           value),
                              M_Fledgling_survival = pull(filter(coucal_vital_rates, 
                                                                 stage == "fledgling" & rate == "survival" & sex == "Male"), 
                                                          value),
                              M_Adult_survival = pull(filter(coucal_vital_rates, 
                                                             stage == "adult" & rate == "survival" & sex == "Male"), 
                                                      value),
            
                              # Define hatching sex ratio
                              HSR = pull(filter(coucal_vital_rates, stage == "HSR"), value),
                              
                              # Define the mating system (h), and clutch size (k)
                              h = pull(filter(coucal_vital_rates, stage == "h"), value),
                              k = pull(filter(coucal_vital_rates, stage == "k"), value))
    
    # Build matrix based on rates specified in the list above
    demographic_matrix <- coucal_matrix(demographic_rates_ = demographic_rates)
    
    # Determine the ASR at the stable stage distribution
    ASR_SSD <- matrix_ASR(M = demographic_matrix,
                          h = demographic_rates$h,
                          HSR = demographic_rates$HSR, 
                          iterations = 100,
                          num_boot = num_boot,
                          species = species,
                          iter_add = 1,
                          n = rep(10, nrow(demographic_matrix)))
    
    # Extract ASR
    ASR_estimate <- ASR_SSD$ASR
    
    # make a list of all the results from this iteration
    bootstrap_results_list <- 
      list(offspring_hazard_lambda = haz_ss_lambda, 
           offspring_hazard_rates = haz_ss_curve,
           coucal_vital_rates = coucal_vital_rates,
           ASR_SSD = ASR_SSD,
           bootstrapped_data = coucal_boot_list)
    
    # extract file directory
    file_directory <- paste0(getwd(), output.directory, "/", bootstrap_name, "_", num_boot, ".Rds") # DEMO "/output/bootstraps/hazard/raw/"
    
    # save the results of the iteration
    save(bootstrap_results_list, file = file_directory)
    
    result = bootstrap_results_list
  
  result
}

hazard_boot_out_wrangle()

If the product of run_bootstrap_juv_hazd_ad_surv_ASR() is saved on disk, this function will tidy up the components of the raw output for subsequent analyses

arguments:

  • species: species name (either “BC” or “WBC”)
  • niter: number of iterations specified in the stochastic simulation
  • output_dir: location in the project directory where the output files are to be found by the function (e.g., “tmp”)
  • rds_file: name of the stochastic simulation rds file saved to disk
# single_model_boot_out_wrangle() 
hazard_boot_out_wrangle <- 
  function(species, niter, output_dir, rds_file){
    
    # specify directory string
    boot_path <- paste0(output_dir, species, rds_file, ".rds")
    
    # load each of the four bootstrap output files
    bootstrap_out <- readRDS(file = boot_path)
    
    # bind all the fledgling daily survival rates output together
    hazard_rates_boot <-
      do.call(rbind, lapply(seq(from = 1, to = niter, by = 1),
                            function(x) bootstrap_out["offspring_hazard_rates", x]$offspring_hazard_rates)) %>%
      arrange(iter, sex, age)
    
    # bind all the ASR estimates output together
    ASR_ests_boot <-
      do.call(rbind, lapply(seq(from = 1, to = niter, by = 1),
                            function(x) bootstrap_out["ASR_SSD", x]$ASR_SSD$ASR)) %>% 
      as.data.frame() %>% 
      mutate(species = species)
    
    # bind all the ASR estimates output together
    vital_rate_ests_boot <-
      do.call(rbind, lapply(seq(from = 1, to = niter, by = 1),
                            function(x) bootstrap_out["coucal_vital_rates", x]$coucal_vital_rates))
    
    # tidy all the bound data together into a mega list of results
    boot_out_tidy <- 
      list(hazard_rates_boot = hazard_rates_boot,
           vital_rate_ests_boot = vital_rate_ests_boot,
           ASR_ests_boot = ASR_ests_boot)
    
    boot_out_tidy
  }

sex_diff_hazard()

takes the list of results from the bootstrap procedure and calculates the sex difference in survival at each life stage for each iteration

arguments:

  • boot_out_list: the object produced by the hazard_boot_out_wrangle() above
  • niter: number of iterations specified in the stochastic simulation
sex_diff_hazard <- function(boot_out_list, niter) {
  
  # make an empty datarame to store the results
  sex_diff_surv_output <- data.frame(Adult = niter,
                                     Fledgling = niter,
                                     Groundling = niter,
                                     Nestling = niter,
                                     ISR = niter,
                                     HSR = niter,
                                     h = niter)
  
  # for loop to go through each iteration and calculate the difference between 
  # female and male survival rates for each stage.
  for(i in 1:niter){
    Adult <- 
      pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "adult" & rate == "survival" & sex == "Female"), value) -
      pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "adult" & rate == "survival" & sex == "Male"), value)
    Fledgling <- 
      pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "fledgling" & rate == "survival" & sex == "Female"), value) -
      pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "fledgling" & rate == "survival" & sex == "Male"), value)
    Groundling <- 
      pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "groundling" & rate == "survival" & sex == "Female"), value) -
      pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "groundling" & rate == "survival" & sex == "Male"), value)
    Nestling <- 
      pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "nestling" & rate == "survival" & sex == "Female"), value) -
      pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "nestling" & rate == "survival" & sex == "Male"), value)
    HSR <- 
      0.5 - pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "HSR"), value)
    h <- 
      1 - pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ], 
                  stage == "h"), value)
    
    sex_diff_surv_output[i, 1] <- ifelse(length(Adult) > 0, Adult, NA)
    sex_diff_surv_output[i, 2] <- ifelse(length(Fledgling) > 0, Fledgling, NA)
    sex_diff_surv_output[i, 3] <- ifelse(length(Groundling) > 0, Groundling, NA)
    sex_diff_surv_output[i, 4] <- ifelse(length(Nestling) > 0, Nestling, NA)
    sex_diff_surv_output[i, 6] <- ifelse(length(HSR) > 0, HSR, NA)
    sex_diff_surv_output[i, 7] <- ifelse(length(h) > 0, h, NA)
    
    
  }
  
  # restructure the output and lable columns
  sex_diff_surv_output <- reshape2::melt(data = sex_diff_surv_output)
  colnames(sex_diff_surv_output) <- c("stage", "difference")
  
  # return the output
  sex_diff_surv_output
}

Wrangling survival data

The following chunks show how we wrangled the raw field data to consolidate the data into the appropriate format for the survival analyses

Offspring

status_dat_all <- 
  # read raw data
  read.csv("data/raw/Coucal_chick_survival_2001-2019_20200129.csv", 
           header = TRUE, stringsAsFactors = FALSE, na.strings = c("", " ", "NA")) %>% 
  
  # rename ring_ID column
  dplyr::rename(ring_ID = Ring_ID) %>% 
  
  # make all entries lower case for consistency
  mutate(Fledged_status = tolower(Fledged.),
         site = tolower(site)) %>% 
  
  # select variables of interest
  dplyr::select(species, ring_ID, lab_no, sex, year, site, nest_ID, pref_age, 
                Fledged_status, postf_age, postf_status, ageC, statusC, 
                lay_date, hatch_order) %>% 
  
  # remove all white space from data
  mutate(across(everything(), ~str_trim(.x))) %>% 
  mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>% 
  
  # specify empty data as NA
  mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>% 
  
  # exclude all individuals that died in the nest
  # filter(Fledged_status == "yes") %>% 
  
  # classify columns
  mutate(sex = as.factor(sex),
         ageC = as.numeric(ageC),
         postf_age = as.numeric(postf_age),
         postf_status = as.numeric(postf_status),
         hatch_order = as.numeric(hatch_order),
         pref_age = as.numeric(pref_age)) %>% 
  
  # remove rows with missing sex, age, and status info
  filter(!is.na(sex) & !is.na(ageC) & !is.na(statusC)) %>% 
  
  # make a unique id for each individual
  mutate(ind_ID = paste(nest_ID, lab_no, ring_ID, sep = "_"),
         
         # create the age of entry into the data (all at age 15)
         entry = 0,
         
         # specify the age of death or censoring
         exit = ageC,
         
         # make the event numeric and specify if 
         # the individual died (1) or was censored (0)
         event = as.numeric(statusC)) %>% 
  
  # consolidate to variables of interest
  dplyr::select(species, ind_ID, nest_ID, year, sex, entry, exit, event, hatch_order)

Adult annual survival

detect_dat <-
  read_xlsx("data/raw/All_coucal_waypoints_2001_2019_20200202.xlsx", na = "NA", col_types = "text") %>% 
  dplyr::select(species, ring_ID, sex, year, site, age_status, date_dec) %>% 
  mutate(month = str_sub(date_dec, start = 5, end = 6),
         day = str_sub(date_dec, start = 7, end = 8)) %>% 
  mutate(date = as.Date(paste(year, month, day, sep = "-"), format = "%Y-%m-%d")) %>% 
  dplyr::select(-month, -day, -date_dec) %>% 
  mutate(across(c("sex", "site", "age_status"), tolower)) %>%
  mutate(sex = ifelse(sex == "female", "F", ifelse(sex == "male", "M", "XXX")),
         age_status = ifelse(age_status == "adult", "A", ifelse(age_status == "juvenile", "J", "XXX")),
         ring_ID = str_replace_all(string = ring_ID, fixed(" "), "")) %>%
  mutate(across(c("species", "ring_ID", "sex", "site", "age_status"), as.factor))

# assess the first and last year of study for each species
detect_dat %>% 
  dplyr::group_by(species) %>% 
  dplyr::summarise(first_year = min(year),
            last_year = max(year))
# A tibble: 3 × 3
  species first_year last_year
  <fct>   <chr>      <chr>    
1 BC      2001       2018     
2 CTC     2016       2019     
3 WBC     2005       2019     
# make annual capture histories for adults
BC_detect_dat_A <-
  detect_dat %>% 
  filter(age_status == "A" & species == "BC") %>% 
  mutate(year = as.numeric(year))

BC_detect_dat_A %>% 
  group_by(ring_ID) %>% 
  dplyr::summarise(n_years = n_distinct(year)) %>% 
  arrange(desc(n_years))
# A tibble: 280 × 2
   ring_ID           n_years
   <fct>               <int>
 1 GN41752                 2
 2 GN52308                 2
 3 GN55932                 2
 4 GN76761                 2
 5 GN76762                 2
 6 bar_tailed_female       1
 7 BCF#34                  1
 8 BCF#35                  1
 9 BCF#37                  1
10 BCF#39                  1
# ℹ 270 more rows
# use the BaSTA function "CensusToCaptHist()" function to convert long format
# encounter histories of each individual, to wide format with 1's and 0's for 
# each year of encounter
BC_detect_dat_A_ch <- 
  CensusToCaptHist(ID = BC_detect_dat_A$ring_ID,
                   d = BC_detect_dat_A$year) %>% 
  dplyr::rename(ring_ID = ID) %>%
  mutate(ID = rownames(.)) %>% 
  left_join(., dplyr::select(BC_detect_dat_A, ring_ID, sex, year), by = "ring_ID") %>% 
  distinct()

Black_Coucal_adult_CJS_ch <- 
  data.frame(ch = apply(BC_detect_dat_A_ch[, 2:19 ] , 1, paste, collapse = "")) %>% 
  bind_cols(., dplyr::select(BC_detect_dat_A_ch, sex, ring_ID)) %>% 
  mutate(across(everything(), ~str_trim(.x))) %>% 
  mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>% 
  mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>% 
  distinct()

WBC_detect_dat_A <-
  detect_dat %>% 
  filter(age_status == "A" & species == "WBC") %>% 
  mutate(year = as.numeric(year))

WBC_detect_dat_A_ch <- 
  CensusToCaptHist(ID = WBC_detect_dat_A$ring_ID,
                   d = WBC_detect_dat_A$year) %>% 
  dplyr::rename(ring_ID = ID) %>%
  mutate(ID = rownames(.)) %>% 
  left_join(., dplyr::select(WBC_detect_dat_A, ring_ID, sex, year), by = "ring_ID") %>% 
  distinct()

White_browed_Coucal_adult_CJS_ch <- 
  data.frame(ch = apply(WBC_detect_dat_A_ch[, 2:16 ] , 1, paste, collapse = "")) %>% 
  bind_cols(., dplyr::select(WBC_detect_dat_A_ch, sex, ring_ID)) %>% 
  mutate(across(everything(), ~str_trim(.x))) %>% 
  mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>% 
  mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>% 
  distinct()

Bootstrapped mark-recapture analysis of White-browed Coucal adults

Here we resample the annual encounter history of adult White-browed Coucals and conduct a CJS mark-recapture analysis to derive the mean estimate and variance in annual apparent survival.

Set the ‘niter’ to run the bootstrap over several iterations (e.g., 10000). For reproducibility, set the seed with the set.seed() function. Set the ‘output.directory’ so that the subsequent chunks can refer to the appropriate folder within the project directory.

niter = 10000
set.seed(14)
output.directory = "output/bootstraps/"

run bootstrap procedure on White-browed Coucals

WBC_adult_survival_CJS_bootstrap_result <-
  pbsapply(1:niter, run_bootstrap_adult_survival_CJS, 
           adult = White_browed_Coucal_adult_CJS_ch,
           species = "WBC",
           iter_add = 1,
           first_year = 2005,
           bootstrap_name = "WBC_boot_one",
           prefix_number = "boot_one_",
           out_directory = "output/bootstraps/")
# save the output
saveRDS(WBC_adult_survival_CJS_bootstrap_result,
        file = "output/bootstraps/WBC_adult_survival_CJS_bootstrap_result.rds")

# tidy it up
WBC_adult_survival_CJS_bootstrap_tidy <- 
  surv_boot_out_wrangle(species = "WBC", niter = niter, output_directory = output.directory)

# save the tidy version
saveRDS(WBC_adult_survival_CJS_bootstrap_tidy,
        file = "output/bootstraps/WBC_adult_survival_CJS_bootstrap_tidy.rds")

Derive the Black Coucal sex-specific survival rates by transforming White-browed Coucal sex-specific survival rates to reflect the ~0.7 ASR that is observed annually

# bootstrapped WBC survival rates
WBC_adult_surival_boot_out <-
  WBC_boot_out$adult_reals_survival_rates_boot %>% 
  dplyr::rename(value = estimate)

# boostrapped BC survival rates (i.e., derived from the WBC rates)
trans_WBC_adult_surival_boot_out <-
  WBC_boot_out$adult_reals_survival_rates_boot %>% 
  dplyr::rename(value = estimate) %>% 
  mutate(value = ifelse(sex == "Female", value * 0.6, value * 1.4))

Parameterizing demographic rates

Each of the following rates has a mean and standard deviation that are used in the stochastic simulation of adult sex ratio.

Mating system

mating_dat <- 
  # read raw data
  read.csv("data/raw/Coucal_Nr_mates_2001_2019_20200123.csv", 
           header = TRUE, stringsAsFactors = FALSE, na.strings = c("", " ", "NA")) %>% 

  # make all entries lower case for consistency
  mutate(site = tolower(site),
         age_status = tolower(age_status),
         sex = tolower(sex)) %>% 

  # specify empty data as NA
  mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>% 
  
  # classify columns
  mutate(sex = as.factor(sex),
         species = as.factor(species),
         Nr_partners = as.numeric(Nr_partners),
         age_status = as.factor(age_status),
         site = as.factor(site)) %>% 
  
  # remove rows with missing sex and age, and post-fledged status info
  filter(!is.na(sex) & !is.na(Nr_partners) & !is.na(species) & species != "CTC") %>% 
  
  # make a unique id for each individual
  mutate(ind_ID = str_replace_all(ring_ID, fixed(" "), "")) %>% 
  mutate(ID_partner.s. = str_replace_all(ID_partner.s., fixed(" "), "")) %>% 

  # make a unique id for each individual
  mutate(ind_ID = str_replace_all(ind_ID, fixed("_"), "")) %>% 
  mutate(ID_partner.s. = str_replace_all(ID_partner.s., fixed("_"), "")) %>% 
  
  
  # consolidate to variables of interest
  dplyr::select(year, species, ind_ID, site, sex, age_status, Nr_partners, ID_partner.s.)

# sex_specific_mating_system <- 
#   mating_dat %>% 
#   group_by(species, sex) %>% 
#   dplyr::summarise(mean_annual_no_mates = mean(Nr_partners, na.rm = TRUE),
#                    var_annual_no_mates = var(Nr_partners, na.rm = TRUE),
#                    median_annual_no_mates = median(Nr_partners, na.rm = TRUE),
#                    sd_annual_no_mates = sd(Nr_partners, na.rm = TRUE),
#                    n = n_distinct(ind_ID), .groups = "drop",
#                    n_years = n_distinct(year)) %>% 
#   mutate(sample_size = paste("n = ", n, sep = ""),
#          species_plot = ifelse(species == "WBC", 1.9, 1.1),
#          species_lab = ifelse(species == "WBC", 1, 2))

sex_specific_mating_system <- 
  mating_dat %>% 
  group_by(ind_ID, sex, species) %>% 
  dplyr::summarise(mean_annual_no_mates_ind = mean(Nr_partners, na.rm = TRUE),
                   n_years = n_distinct(year)) %>% 
  mutate(mu_i = ifelse(mean_annual_no_mates_ind <= 1, 1, mean_annual_no_mates_ind)) %>% 
  group_by(species, sex) %>% 
  dplyr::summarise(mean_annual_no_mates = mean(mu_i, na.rm = TRUE),
                   var_annual_no_mates = var(mu_i, na.rm = TRUE),
                   median_annual_no_mates = median(mu_i, na.rm = TRUE),
                   sd_annual_no_mates = sd(mu_i, na.rm = TRUE),
                   n = n_distinct(ind_ID), .groups = "drop") %>% 
  mutate(sample_size = paste("n = ", n, sep = ""),
         species_plot = ifelse(species == "WBC", 1.9, 1.1),
         species_lab = ifelse(species == "WBC", 1, 2))

# To obtain a female-based h index for each population the inverse of the mean 
# mu is calculated (i.e., Eq. 2 in manuscript)
BC_h <- 
  as.numeric(sex_specific_mating_system[which(
    sex_specific_mating_system$sex == "female" &
    sex_specific_mating_system$species == "BC"), "mean_annual_no_mates"])
WBC_h <- 
  as.numeric(sex_specific_mating_system[which(
    sex_specific_mating_system$sex == "female" &
    sex_specific_mating_system$species == "WBC"), "mean_annual_no_mates"])

BC_h_sd <- 
  as.numeric(sex_specific_mating_system[which(
    sex_specific_mating_system$sex == "female" &
      sex_specific_mating_system$species == "BC"), "sd_annual_no_mates"])

WBC_h_sd <- 
  as.numeric(sex_specific_mating_system[which(
    sex_specific_mating_system$sex == "female" &
      sex_specific_mating_system$species == "WBC"), "sd_annual_no_mates"])

coucal_mating_system <- 
  data.frame(trait = c("mating_system"),
             species = c("BC", "WBC"),
             mean = c(BC_h, WBC_h),
             sd = c(BC_h_sd, WBC_h_sd))

mating_dat$species <- 
  factor(mating_dat$species ,
         levels = c("BC",
                    "WBC"))

mating_dat_plotting <-
  mating_dat %>% 
  mutate(species_dummy = ifelse(species == "WBC", 1, 0)) %>%
  mutate(species_plot = ifelse(species == "WBC", 2.1, 0.9))

mating_system_plot <-
  ggplot2::ggplot() +
  theme_bw() +
  geom_jitter(aes(y = Nr_partners, x = species_plot, 
                  color = sex),
              data = mating_dat_plotting, 
              width = 0.1, height = 0.1,
              alpha = 0.75, 
              fill = "white", shape = 16) +
  geom_pointrange(data = sex_specific_mating_system, 
    aes(y = mean_annual_no_mates, x = species_plot, 
        ymin = (mean_annual_no_mates-sd_annual_no_mates), 
        ymax = (mean_annual_no_mates+sd_annual_no_mates),
        fill = sex),
    size = 0.8, shape = 21, color = "black") +
  geom_text(aes(y = 0.2, x = species_lab, label = sample_size), 
            data = sex_specific_mating_system,
            size = 3, family = "Verdana") +
  geom_hline(yintercept = 1.5, linetype = "dashed",
             alpha = 0.5, color = "grey20") +
  facet_grid(. ~ sex, labeller = as_labeller(sex_names)) +
  luke_theme +
  theme(legend.position = "none",
        axis.title.x = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size = 12, face = "italic")) +
  scale_x_continuous(limits = c(0.5, 2.5), 
                   breaks = c(1, 2), labels = species_names) +
  scale_y_continuous(limits = c(0, 5.5), expand = c(0, 0), 
                     breaks = c(1, 2, 3, 4, 5),
                     labels = c(1, 2, 3, 4, 5)) +
  ylab("Number of unique mates\nper year (± 1 SD)") +
  scale_color_manual(values = plot_palette_sex) +
  scale_fill_manual(values = plot_palette_sex) +
  annotate(geom = "text", y = 1.25, x = 1.5,
           label = "monogamy", family = "Verdana", 
           color = "black", size = 3, fontface = 'italic', hjust = 0.5) +
  annotate(geom = "text", y = 1.75, x = 1.5,
           label = "polygamy", family = "Verdana",
           color = "black", size = 3, fontface = 'italic', hjust = 0.5)

Clutch size

# import raw csv data into R
egg_data <-
  read_xls("data/raw/Egg_surv_data_2001_2019_20210524.xls", na = "NA", col_types = "text") %>% 
  dplyr::select(species, nest_ID, year, `Total clutch size`, all_chicks, 
                nst_suc_sho, nest_success, `male_chi+eggs`, `female_chi+eggs`,
                `unkn_sex_eggs+chiks`) %>% 
  filter(species != "CTC") %>% 
  filter(!is.na(nest_success)) %>% 
  dplyr::rename(n_eggs = `Total clutch size`,
                n_chicks = all_chicks,
                male_eggs = `male_chi+eggs`, 
                female_eggs = `female_chi+eggs`,
                unknown_eggs = `unkn_sex_eggs+chiks`) %>% 
  dplyr::mutate(n_eggs = na_if(n_eggs, "?"),
                n_chicks = na_if(n_chicks, "?")) %>% 
  dplyr::mutate(n_eggs = as.numeric(n_eggs),
                n_chicks = as.numeric(n_chicks),
                male_eggs = as.numeric(male_eggs),
                female_eggs = as.numeric(female_eggs),
                unknown_eggs = as.numeric(unknown_eggs)) %>% 
  filter(n_eggs >= n_chicks)

#### clutch size summary ----
coucal_clutch_size <- 
  data.frame(trait = c("clutch_size"),
             species = c("BC", "WBC"),
             mean = c(filter(egg_data, species == "BC") %>% 
                                    summarise(mean_n_eggs = mean(n_eggs, na.rm = TRUE)) %>% 
                                    pull(mean_n_eggs),
                                  filter(egg_data, species == "WBC") %>% 
                                    summarise(mean_n_eggs = mean(n_eggs, na.rm = TRUE)) %>% 
                                    pull(mean_n_eggs)),
             sd = c(filter(egg_data, species == "BC") %>% 
                                  summarise(mean_n_eggs = sd(n_eggs, na.rm = TRUE)) %>% 
                                  pull(mean_n_eggs),
                                filter(egg_data, species == "WBC") %>% 
                                  summarise(mean_n_eggs = sd(n_eggs, na.rm = TRUE)) %>% 
                                  pull(mean_n_eggs)),
             n_nests = c(filter(egg_data, species == "BC") %>% 
                           summarise(n_nests = n_distinct(nest_ID)) %>% 
                           pull(n_nests),
                         filter(egg_data, species == "WBC") %>% 
                           summarise(n_nests = n_distinct(nest_ID)) %>% 
                           pull(n_nests)),
             n_years = c(filter(egg_data, species == "BC") %>% 
                           summarise(n_nests = n_distinct(year)) %>% 
                           pull(n_nests),
                         filter(egg_data, species == "WBC") %>% 
                           summarise(n_nests = n_distinct(year)) %>% 
                           pull(n_nests)))

Egg hatchability

# import raw csv data into R
egg_data <-
  read_xls("data/raw/Egg_surv_data_2001_2019_20210524.xls", 
           na = "NA", col_types = "text") %>% 
  dplyr::select(species, nest_ID, year, `Total clutch size`, all_chicks, 
                nst_suc_sho, nest_success, `male_chi+eggs`, `female_chi+eggs`,
                `unkn_sex_eggs+chiks`) %>% 
  filter(species != "CTC") %>% 
  filter(!is.na(nest_success)) %>% 
  dplyr::rename(n_eggs = `Total clutch size`,
                n_chicks = all_chicks,
                male_eggs = `male_chi+eggs`, 
                female_eggs = `female_chi+eggs`,
                unknown_eggs = `unkn_sex_eggs+chiks`) %>% 
  dplyr::mutate(n_eggs = na_if(n_eggs, "?"),
                n_chicks = na_if(n_chicks, "?")) %>%
  dplyr::mutate(n_eggs = as.numeric(n_eggs),
                n_chicks = as.numeric(n_chicks),
                male_eggs = as.numeric(male_eggs),
                female_eggs = as.numeric(female_eggs),
                unknown_eggs = as.numeric(unknown_eggs)) %>% 
  filter(n_eggs >= n_chicks) %>% 
  arrange(nest_ID)

#### egg hatchability summary ----
egg_surv_data <-
  egg_data %>% 
  dplyr::mutate(n_dead_eggs = n_eggs - n_chicks,
                n_alive_eggs = n_chicks) %>% 
  arrange(nest_ID)

#### sample size summary ----
egg_surv_data %>% 
  dplyr::group_by(species) %>% 
  dplyr::summarise(n_nests = n_distinct(nest_ID))
# A tibble: 2 × 2
  species n_nests
  <chr>     <int>
1 BC          271
2 WBC         213
BC_egg_survival <- 
  lme4::glmer(cbind(n_alive_eggs, n_dead_eggs) ~ 
                1 + 
                (1| nest_ID), 
              data = filter(egg_surv_data, species == "BC"), 
              family = binomial)

WBC_egg_survival <- 
  lme4::glmer(cbind(n_alive_eggs, n_dead_eggs) ~ 
                1 + 
                (1| nest_ID), 
              data = filter(egg_surv_data, species == "WBC"), 
              family = binomial)

# # run tidy bootstrap to obtain model diagnostics
tidy_BC_egg_survival <-
  tidy(BC_egg_survival, conf.int = TRUE, conf.method = "boot", nsim = 1000) # DEMO 1000

tidy_WBC_egg_survival <-
  tidy(WBC_egg_survival, conf.int = TRUE, conf.method = "boot", nsim = 1000) # DEMO 1000

coucal_egg_survival <- 
  data.frame(trait = c("egg_survival"),
             species = c("BC", "WBC"),
             mean = c(invlogit(model_parameters(BC_egg_survival)$Coefficient)[1],
                      invlogit(model_parameters(WBC_egg_survival)$Coefficient)[1]),
             CI_low = c(invlogit(model_parameters(BC_egg_survival)$CI_low)[1],
                        invlogit(model_parameters(WBC_egg_survival)$CI_low)[1]),
             CI_high = c(invlogit(model_parameters(BC_egg_survival)$CI_high)[1],
                         invlogit(model_parameters(WBC_egg_survival)$CI_high)[1]),
             n_nests = c(filter(egg_data, species == "BC") %>% 
                           summarise(n_nests = n_distinct(nest_ID)) %>% 
                           pull(n_nests),
                         filter(egg_data, species == "WBC") %>% 
                           summarise(n_nests = n_distinct(nest_ID)) %>% 
                           pull(n_nests)),
             n_years = c(filter(egg_data, species == "BC") %>% 
                           summarise(n_nests = n_distinct(year)) %>% 
                           pull(n_nests),
                         filter(egg_data, species == "WBC") %>% 
                           summarise(n_nests = n_distinct(year)) %>% 
                           pull(n_nests))) %>% 
  mutate(sd = ifelse(!is.na(CI_low), 
                         approx_sd(x1 = CI_low, x2 = CI_high),
                         CI_low))

Hatching sex ratio

egg_data <-
  read_xls("data/raw/Egg_surv_data_2001_2019_20210524.xls", 
           na = "NA", col_types = "text") %>% 
  dplyr::select(species, nest_ID, year, `Total clutch size`, all_chicks, 
                nst_suc_sho, nest_success, `male_chi+eggs`, `female_chi+eggs`,
                `unkn_sex_eggs+chiks`) %>% 
  filter(species != "CTC") %>% 
  filter(!is.na(nest_success)) %>% 
  dplyr::rename(n_eggs = `Total clutch size`,
                n_chicks = all_chicks,
                male_eggs = `male_chi+eggs`, 
                female_eggs = `female_chi+eggs`,
                unknown_eggs = `unkn_sex_eggs+chiks`) %>% 
  dplyr::mutate(n_eggs = na_if(n_eggs, "?"),
                n_chicks = na_if(n_chicks, "?")) %>% 
  dplyr::mutate(n_eggs = as.numeric(n_eggs),
                n_chicks = as.numeric(n_chicks),
                male_eggs = as.numeric(male_eggs),
                female_eggs = as.numeric(female_eggs),
                unknown_eggs = as.numeric(unknown_eggs)) %>% 
  filter(n_eggs >= n_chicks)

#### hatching sex ratio summary ----
egg_HSR_data <- 
  egg_data %>% 
  dplyr::filter(n_eggs == (male_eggs + female_eggs + unknown_eggs)) %>% 
  dplyr::filter(unknown_eggs == 0) %>% 
  mutate(HSR = male_eggs/(male_eggs + female_eggs + unknown_eggs)) %>%
  mutate(species_plot = ifelse(species == "WBC", 2.2, 0.8),
         species_plot2 = ifelse(species == "WBC", 2, 1))

egg_HSR_data %>% 
  dplyr::group_by(species) %>% 
  dplyr::summarise(n_distinct(nest_ID))
# A tibble: 2 × 2
  species `n_distinct(nest_ID)`
  <chr>                   <int>
1 BC                         70
2 WBC                        51
mod_HSR_BC <- 
  lme4::glmer(cbind(male_eggs, female_eggs) ~ 
                1 + 
                (1| nest_ID), 
              data = filter(egg_HSR_data, species == "BC"), 
              family = binomial)

mod_HSR_WBC <- 
  lme4::glmer(cbind(male_eggs, female_eggs) ~ 
                1 + 
                (1| nest_ID), 
              data = filter(egg_HSR_data, species == "WBC"), 
              family = binomial)

# # run tidy bootstrap to obtain model diagnostics
tidy_mod_HSR_BC <-
  tidy(mod_HSR_BC, conf.int = TRUE, conf.method = "boot", nsim = 1000) # DEMO 1000

tidy_mod_HSR_WBC <-
  tidy(mod_HSR_WBC, conf.int = TRUE, conf.method = "boot", nsim = 1000) # DEMO 1000

coucal_HSR <- 
  data.frame(trait = "hatching_sex_ratio",
             species = c("BC", "WBC"),
             mean = c(invlogit(model_parameters(mod_HSR_BC)$Coefficient)[1],
                          invlogit(model_parameters(mod_HSR_WBC)$Coefficient)[1]),
             CI_low = c(invlogit(model_parameters(mod_HSR_BC)$CI_low)[1],
                        invlogit(model_parameters(mod_HSR_WBC)$CI_low)[1]),
             CI_high = c(invlogit(model_parameters(mod_HSR_BC)$CI_high)[1],
                         invlogit(model_parameters(mod_HSR_WBC)$CI_high)[1]),
             n_nests = c(filter(egg_HSR_data, species == "BC") %>% 
                           summarise(n_nests = n_distinct(nest_ID)) %>% 
                           pull(n_nests),
                         filter(egg_HSR_data, species == "WBC") %>% 
                           summarise(n_nests = n_distinct(nest_ID)) %>% 
                           pull(n_nests)),
             n_years = c(filter(egg_HSR_data, species == "BC") %>% 
                           summarise(n_nests = n_distinct(year)) %>% 
                           pull(n_nests),
                         filter(egg_HSR_data, species == "WBC") %>% 
                           summarise(n_nests = n_distinct(year)) %>% 
                           pull(n_nests))) %>% 
  mutate(sd = ifelse(!is.na(CI_low), 
                         approx_sd(x1 = CI_low, x2 = CI_high),
                         CI_low),
         species_plot = ifelse(species == "WBC", 1.8, 1.2))

HSR_plot <-
  ggplot2::ggplot() +
  geom_boxplot(data = egg_HSR_data,
               aes(x = species_plot, y = HSR,
                   group = species, fill = species),
               color = "grey50",
               width = 0.05, alpha = 0.5,
               position = position_dodge(width = 0)) +
  geom_errorbar(data = coucal_HSR,
                aes(x = species_plot, ymax = CI_high, ymin = CI_low),
                alpha = 1, color = "black", width = 0.05, lwd = 0.5) +
  geom_point(data = coucal_HSR,
             aes(x = species_plot, y = mean, fill = species),
             lwd = 3, shape = 21, color= "black") +
  geom_jitter(data = egg_HSR_data,
              aes(x = species_plot2, y = HSR,
                  group = species,
                  fill = species, color = species),
              width = 0.02, alpha = 0.2, shape = 19) +
  luke_theme +
  theme(legend.position = "none",
        panel.border = element_blank(),
        strip.background = element_blank(),
        axis.title.x = element_blank(),
        axis.text.x = element_text(12),

        legend.background = element_blank(),
        panel.grid.major.y = element_line(colour = "grey70", size = 0.1),
        plot.title = element_text(hjust = 0.5, face = "italic")) +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.25)) +
  scale_x_discrete(labels = c("1" = "Black Coucal",
                              "2" = "White-browed Coucal")) +
  ylab(expression(paste("Hatching sex ratio (prop. male)" %+-%  "95% CI", sep = ""))) +
  xlab("Species") +
  scale_color_manual(values = c(brewer.pal(6, "Dark2")[1], brewer.pal(6, "Set1")[1])) +
  scale_fill_manual(values = c(brewer.pal(6, "Dark2")[1], brewer.pal(6, "Set1")[1])) +
  ggtitle("Black Coucal           White-browed\n                                Coucal")

Age at fledging (transition from nestling to groundling)

fledge_dat <- 
  read.csv("data/raw/Coucal_chick_survival_2001-2019_20200129.csv", 
           header = TRUE, stringsAsFactors = FALSE, na.strings = c("", " ", "NA")) %>%
  dplyr::select(species, sex, Fledge_age, nest_ID, year, lab_no, Ring_ID, site) %>% 
  filter(site == "Kapunga" & species != "CTC") %>% 
  mutate(Nst_No = str_replace_all(string = nest_ID, fixed(" "), ""),
         Ring_ID = str_replace_all(string = Ring_ID, fixed(" "), ""),
         lab_no = str_replace_all(string = lab_no, fixed(" "), "")) %>% 
  filter(!is.na(Fledge_age)) %>% 
  mutate(Ind_ID = paste(Nst_No, Ring_ID, lab_no, sep = "_")) %>%
  mutate(sex_plot = ifelse(sex == "M", 2.2, 0.8))
 
#### Modeling (BC) ----
# "Fledge_age" as dependent variable, interaction with sex and species
mod_fledge_age_BC <- 
  lmer(Fledge_age ~ sex + 
         (1 | nest_ID) + (1 | year), 
       data = filter(fledge_dat, species == "BC"))

new_data_BC <- 
  expand.grid(sex = c("F","M"))

predict_BC <- 
  predict(mod_fledge_age_BC, 
          newdata = new_data_BC, 
          re.form = NA, 
          se.fit = TRUE, 
          nsim = 1000)

#### Modeling (WBC) ----
# "Fledge_age" as dependent variable, interaction with sex and species
mod_fledge_age_WBC <- 
  lmer(Fledge_age ~ sex + 
         (1 | nest_ID) + (1 | year), 
       data = filter(fledge_dat, species == "WBC"))

new_data_WBC <- 
  expand.grid(sex = c("F","M"))

predict_WBC <- 
  predict(mod_fledge_age_WBC, 
          newdata = new_data_WBC, 
          re.form = NA, 
          se.fit = TRUE, 
          nsim = 1000)

# tidy up estimates
coucal_fledge_age <- 
  data.frame(trait = c("fledge_age"),
             species = c("BC", "BC", "WBC", "WBC"),
             sex = c("F", "M", "F", "M"),
             mean = c(predict_BC$fit[1],
                      predict_BC$fit[2],
                      predict_WBC$fit[1],
                      predict_WBC$fit[2]),
             CI_low = c(predict_BC$ci.fit[1, 1],
                        predict_BC$ci.fit[1, 2],
                        predict_WBC$ci.fit[1, 1],
                        predict_WBC$ci.fit[1, 2]),
             CI_high = c(predict_BC$ci.fit[2, 1],
                         predict_BC$ci.fit[2, 2],
                         predict_WBC$ci.fit[2, 1],
                         predict_WBC$ci.fit[2, 2]),
             n_inds = c(filter(fledge_dat, species == "BC" & sex == "F") %>% 
                          summarise(n_ = n_distinct(Ind_ID)) %>% 
                          pull(n_),
                        filter(fledge_dat, species == "BC" & sex == "M") %>% 
                          summarise(n_ = n_distinct(Ind_ID)) %>% 
                          pull(n_),
                        filter(fledge_dat, species == "WBC" & sex == "F") %>% 
                          summarise(n_ = n_distinct(Ind_ID)) %>% 
                          pull(n_),
                        filter(fledge_dat, species == "BC" & sex == "M") %>% 
                          summarise(n_ = n_distinct(Ind_ID)) %>% 
                          pull(n_)),
             n_nests = c(filter(fledge_dat, species == "BC" & sex == "F") %>% 
                           summarise(n_ = n_distinct(Nst_No)) %>% 
                           pull(n_),
                         filter(fledge_dat, species == "BC" & sex == "M") %>% 
                           summarise(n_ = n_distinct(Nst_No)) %>% 
                           pull(n_),
                         filter(fledge_dat, species == "WBC" & sex == "F") %>% 
                           summarise(n_ = n_distinct(Nst_No)) %>% 
                           pull(n_),
                         filter(fledge_dat, species == "WBC" & sex == "M") %>% 
                           summarise(n_ = n_distinct(Nst_No)) %>% 
                           pull(n_)),
             n_years = c(filter(fledge_dat, species == "BC" & sex == "F") %>% 
                           summarise(n_ = n_distinct(year)) %>% 
                           pull(n_),
                         filter(fledge_dat, species == "BC" & sex == "M") %>% 
                           summarise(n_ = n_distinct(year)) %>% 
                           pull(n_),
                         filter(fledge_dat, species == "WBC" & sex == "F") %>% 
                           summarise(n_ = n_distinct(year)) %>% 
                           pull(n_),
                         filter(fledge_dat, species == "WBC" & sex == "M") %>% 
                           summarise(n_ = n_distinct(year)) %>% 
                           pull(n_))) %>% 
  mutate(sd = ifelse(!is.na(CI_low), 
                         approx_sd(x1 = CI_low, x2 = CI_high),
                         CI_low))

Age at first flight (transition from groundling to fledgling)

flight_dat <- 
  read.csv("data/raw/coucals_first_flights_data_per_individual2020.csv", 
           header = TRUE, stringsAsFactors = FALSE, na.strings = c("", " ", "NA")) %>% 
  mutate(Nst_No = str_replace_all(string = Nst_No, fixed(" "), ""),
         Ind_ID = str_replace_all(string = Ind_ID, fixed(" "), "")) %>%
  dplyr::rename(sex = Sex,
                nest_ID = Nst_No,
                year = Year,
                species = Spp,
                flight_age = days_since_fledging,
                Ring_ID = Ind_ID) %>% 
  mutate(sex_plot = ifelse(sex == "Male", 2.2, 0.8))

#### Modeling (WBC) ----
# "flight_age" as dependent variable, 
# "sex", "fledge_age", and "fledge_mass" as independent
mod_flight_age_WBC <-
  lmer(flight_age ~ sex + Fledge_age + Fledge_mass +
         (1 | nest_ID) + (1 | year), 
       data = filter(flight_dat, species == "WBC"))

new_data_WBC <- 
  expand.grid(sex = c("Female","Male"),
              Fledge_age = mean(flight_dat[flight_dat$species == "WBC",]$Fledge_age),
              Fledge_mass = mean(flight_dat[flight_dat$species == "WBC",]$Fledge_mass))

predict_WBC <- 
  predict(mod_flight_age_WBC, 
        newdata = new_data_WBC, 
        re.form = NA, 
        se.fit = TRUE, 
        nsim = 1000)

#### Modeling (BC) ----
# "flight_age" as dependent variable, 
# "sex", "fledge_age", and "fledge_mass" as independent
mod_flight_age_BC <- 
  lmer(flight_age ~ sex + Fledge_age + Fledge_mass +
         (1 | nest_ID) + (1 | year), 
       data = filter(flight_dat, species == "BC"))

new_data_BC <- 
  expand.grid(sex = c("Female","Male"),
              Fledge_age = mean(flight_dat[flight_dat$species == "BC",]$Fledge_age),
              Fledge_mass = mean(flight_dat[flight_dat$species == "BC",]$Fledge_mass))

predict_BC <- 
  predict(mod_flight_age_BC, 
        newdata = new_data_BC, 
        re.form = NA, 
        se.fit = TRUE, 
        nsim = 1000)

coucal_flight_age <- 
  data.frame(trait = c("flight_age"),
             species = c("BC", "BC", "WBC", "WBC"),
             sex = c("F", "M", "F", "M"),
             mean = c(predict_BC$fit[1],
                      predict_BC$fit[2],
                      predict_WBC$fit[1],
                      predict_WBC$fit[2]),
             CI_low = c(predict_BC$ci.fit[1, 1],
                        predict_BC$ci.fit[1, 2],
                        predict_WBC$ci.fit[1, 1],
                        predict_WBC$ci.fit[1, 2]),
             CI_high = c(predict_BC$ci.fit[2, 1],
                        predict_BC$ci.fit[2, 2],
                        predict_WBC$ci.fit[2, 1],
                        predict_WBC$ci.fit[2, 2]),
             n_inds = c(filter(flight_dat, species == "BC" & sex == "Female") %>% 
                          summarise(n_ = n_distinct(Ring_ID)) %>% 
                          pull(n_),
                        filter(flight_dat, species == "BC" & sex == "Male") %>% 
                          summarise(n_ = n_distinct(Ring_ID)) %>% 
                          pull(n_),
                        filter(flight_dat, species == "WBC" & sex == "Female") %>% 
                          summarise(n_ = n_distinct(Ring_ID)) %>% 
                          pull(n_),
                        filter(flight_dat, species == "BC" & sex == "Male") %>% 
                          summarise(n_ = n_distinct(Ring_ID)) %>% 
                          pull(n_)),
             n_nests = c(filter(flight_dat, species == "BC" & sex == "Female") %>% 
                           summarise(n_ = n_distinct(nest_ID)) %>% 
                           pull(n_),
                         filter(flight_dat, species == "BC" & sex == "Male") %>% 
                           summarise(n_ = n_distinct(nest_ID)) %>% 
                           pull(n_),
                         filter(flight_dat, species == "WBC" & sex == "Female") %>% 
                           summarise(n_ = n_distinct(nest_ID)) %>% 
                           pull(n_),
                         filter(flight_dat, species == "WBC" & sex == "Male") %>% 
                           summarise(n_ = n_distinct(nest_ID)) %>% 
                           pull(n_)),
             n_years = c(filter(flight_dat, species == "BC" & sex == "Female") %>% 
                           summarise(n_ = n_distinct(year)) %>% 
                           pull(n_),
                         filter(flight_dat, species == "BC" & sex == "Male") %>% 
                           summarise(n_ = n_distinct(year)) %>% 
                           pull(n_),
                         filter(flight_dat, species == "WBC" & sex == "Female") %>% 
                           summarise(n_ = n_distinct(year)) %>% 
                           pull(n_),
                         filter(flight_dat, species == "WBC" & sex == "Male") %>% 
                           summarise(n_ = n_distinct(year)) %>% 
                           pull(n_))) %>% 
  mutate(sd = ifelse(!is.na(CI_low), 
                     approx_sd(x1 = CI_low, x2 = CI_high),
                     CI_low))

Bind all parameters together

parameter_distributions <- 
  bind_rows(coucal_egg_survival,
            coucal_clutch_size,
            coucal_HSR,
            coucal_mating_system,
            coucal_fledge_age,
            coucal_flight_age) %>% 
  dplyr::select(trait, species, sex, mean, sd)

Stochastic Demographic Simulation

niter = 10000
set.seed(14)

Black Coucals

pull species- and sex-specific means and sds for each parameter

k_dist_BC = c(pull(filter(coucal_clutch_size, species == "BC"), mean), 
              pull(filter(coucal_clutch_size, species == "BC"), sd))

HSR_dist_BC = c(pull(filter(coucal_HSR, species == "BC"), mean), 
                pull(filter(coucal_HSR, species == "BC"), sd))

h_dist_BC = c(pull(filter(coucal_mating_system, species == "BC"), mean), 
              pull(filter(coucal_mating_system, species == "BC"), sd))

egg_surv_dist_BC = c(pull(filter(coucal_egg_survival, species == "BC"), mean),
                     pull(filter(coucal_egg_survival, species == "BC"), sd))

fledge_age_distF_BC = c(pull(filter(coucal_fledge_age, species == "BC" & sex == "F"), 
                             mean), 
                        pull(filter(coucal_fledge_age, species == "BC" & sex == "F"), 
                             sd))

fledge_age_distM_BC = c(pull(filter(coucal_fledge_age, species == "BC" & sex == "M"), 
                             mean), 
                        pull(filter(coucal_fledge_age, species == "BC" & sex == "M"), 
                             sd))

flight_age_distF_BC = c(pull(filter(coucal_flight_age, species == "BC" & sex == "F"), 
                             mean), 
                        pull(filter(coucal_flight_age, species == "BC" & sex == "F"), 
                             sd))

flight_age_distM_BC = c(pull(filter(coucal_flight_age, species == "BC" & sex == "M"), 
                             mean), 
                        pull(filter(coucal_flight_age, species == "BC" & sex == "M"), 
                             sd))

run stochastic demographic simulation

BC_hazard_ASR_bootstrap_result_w_trans_WBC_ad_surv_stoc <-
  pbsapply(1:niter, run_bootstrap_juv_hazd_ad_surv_ASR,
           offspring = status_dat_all %>% filter(species == "BC"),
           adult_surival_boot_out = trans_WBC_adult_surival_boot_out,
           bootstrap_name = "BC_boot_w_trans_WBC_ad_surv_stoc",
           species = "BC",
           iter_add = 1,
           prefix_number = "boot_w_trans_WBC_ad_surv_stoc_no",
           max_time = 70,
           niter = niter,
           alpha_value = 1.4,
           k_dist = k_dist_BC,
           HSR_dist = HSR_dist_BC,
           h_dist = h_dist_BC,
           egg_surv_dist = egg_surv_dist_BC,
           fledge_age_distF = fledge_age_distF_BC,
           fledge_age_distM = fledge_age_distM_BC,
           flight_age_distF = flight_age_distF_BC,
           flight_age_distM = flight_age_distM_BC, 
           output.directory = "/output/bootstraps/hazard/raw")

White-browed Coucals

pull species- and sex-specific means and sds for each parameter

k_dist_WBC = c(pull(filter(coucal_clutch_size, species == "WBC"), mean), 
               pull(filter(coucal_clutch_size, species == "WBC"), sd))

HSR_dist_WBC = c(pull(filter(coucal_HSR, species == "WBC"), mean), 
                 pull(filter(coucal_HSR, species == "WBC"), sd))

h_dist_WBC = c(pull(filter(coucal_mating_system, species == "WBC"), mean), 
               pull(filter(coucal_mating_system, species == "WBC"), sd))

egg_surv_dist_WBC = c(pull(filter(coucal_egg_survival, species == "WBC"), mean),
                      pull(filter(coucal_egg_survival, species == "WBC"), sd))

fledge_age_distF_WBC = c(pull(filter(coucal_fledge_age, species == "WBC" & sex == "F"), 
                              mean), 
                         pull(filter(coucal_fledge_age, species == "WBC" & sex == "F"), 
                              sd))

fledge_age_distM_WBC = c(pull(filter(coucal_fledge_age, species == "WBC" & sex == "M"), 
                              mean), 
                         pull(filter(coucal_fledge_age, species == "WBC" & sex == "M"), 
                              sd))

flight_age_distF_WBC = c(pull(filter(coucal_flight_age, species == "WBC" & sex == "F"), 
                              mean), 
                         pull(filter(coucal_flight_age, species == "WBC" & sex == "F"), 
                              sd))

flight_age_distM_WBC = c(pull(filter(coucal_flight_age, species == "WBC" & sex == "M"),
                              mean), 
                         pull(filter(coucal_flight_age, species == "WBC" & sex == "M"), 
                              sd))

run stochastic demographic simulation

WBC_hazard_ASR_bootstrap_result_w_WBC_ad_surv_stoc <-
  pbsapply(1:niter, run_bootstrap_juv_hazd_ad_surv_ASR,
           offspring = status_dat_all %>% filter(species == "WBC"), 
           adult_surival_boot_out = WBC_adult_surival_boot_out,
           bootstrap_name = "WBC_boot_w_WBC_ad_surv_stoc",
           species = "WBC",
           iter_add = 1,
           prefix_number = "boot_w_WBC_ad_surv_stoc", 
           max_time = 70,
           niter = niter,
           alpha_value = 1.4,
           k_dist = k_dist_WBC,
           HSR_dist = HSR_dist_WBC,
           h_dist = h_dist_WBC,
           egg_surv_dist = egg_surv_dist_WBC,
           fledge_age_distF = fledge_age_distF_WBC,
           fledge_age_distM = fledge_age_distM_WBC,
           flight_age_distF = flight_age_distF_WBC,
           flight_age_distM = flight_age_distM_WBC,
           output.directory = "/output/bootstraps/hazard/raw")

Save simulation output

saveRDS(object = BC_hazard_ASR_bootstrap_result_w_trans_WBC_ad_surv_stoc, 
        file = "output/bootstraps/hazard/cooked/BC_hazard_ASR_bootstrap_result_w_trans_WBC_ad_surv_stoc.rds")

saveRDS(object = WBC_hazard_ASR_bootstrap_result_w_WBC_ad_surv_stoc, 
        file = "output/bootstraps/hazard/cooked/WBC_hazard_ASR_bootstrap_result_w_WBC_ad_surv_stoc.rds")

clean up the output from the bootstrap procedure and save as rds

BC_hazard_rate_boot_tidy_trans_no_imm <- 
  hazard_boot_out_wrangle(species = "BC", niter = 1000, 
                          output_dir = "output/bootstraps/hazard/cooked/",
                          rds_file = 
                           "_hazard_ASR_bootstrap_result_w_trans_WBC_ad_surv_stoc_no_imm")

WBC_hazard_rate_boot_tidy_no_imm <- 
  hazard_boot_out_wrangle(species = "WBC", niter = 1000, 
                          output_dir = "output/bootstraps/hazard/cooked/",
                          rds_file = 
                            "_hazard_ASR_bootstrap_result_w_WBC_ad_surv_stoc_no_imm")

save model output

saveRDS(object = BC_hazard_rate_boot_tidy_trans_no_imm, 
        file = 
          "output/bootstraps/hazard/cooked/BC_haz_sur_ASR_boot_tidy_stoc_trans_no_imm.rds")

saveRDS(object = WBC_hazard_rate_boot_tidy_no_imm, 
        file = 
          "output/bootstraps/hazard/cooked/WBC_haz_sur_ASR_boot_tidy_stoc_no_imm.rds")

Sex-biased Demographic stages

# calculate the sex differences in stage specific rates
BC_sex_diff_hazard_output <- 
  sex_diff_hazard(boot_out_list = BC_hazard_rate_boot_tidy, 
                  niter = 1000) %>% 
  mutate(species = "BC")

WBC_sex_diff_hazard_output <- 
  sex_diff_hazard(boot_out_list = WBC_hazard_rate_boot_tidy, 
                  niter = 1000) %>% 
  mutate(species = "WBC")

# consolidate results
sex_diff_survival_output <- 
  dplyr::bind_rows(BC_sex_diff_hazard_output,
                   WBC_sex_diff_hazard_output) %>% 
  dplyr::filter(stage %in% c("HSR", "Nestling", "Groundling", 
                             "Fledgling", "Adult")) %>% 
  mutate(stage = factor(stage, 
                        levels = c("HSR", "Nestling", "Groundling", 
                                   "Fledgling", "Adult")),
         species = factor(species, 
                          levels = c("BC", "WBC")),
         difference = -difference)

# check that there are no nosensical values
sex_diff_survival_output %>% 
  filter(difference < -1 & difference > 1)

ggplot(data = filter(sex_diff_survival_output, stage != "h")) +
  geom_histogram(aes(difference), binwidth = 0.01) +
  facet_grid(stage ~ species) +
  scale_x_continuous(limits = c(-0.5, 0.5)) +
  geom_vline(xintercept = 0)

# calculate some summary statistics
sex_diff_survival_summary <- 
  sex_diff_survival_output %>%
  dplyr::group_by(stage, species) %>%
  dplyr::summarise(avg = mean(difference, na.rm = TRUE),
                   median = median(difference, na.rm = TRUE),
                   var = var(difference, na.rm = TRUE),
                   max = max(difference, na.rm = TRUE),
                   min = min(difference, na.rm = TRUE),
                   sd = sd(difference, na.rm = TRUE))

# specify custom color palette to distingush first-year stages 
# (i.e. chicks and juveniles) from adults
cbPalette <- c("#D9D9D9", "#D9D9D9", "#D9D9D9", "#D9D9D9", "#A6A6A6")

species_names <- 
  c('BC' = "Black Coucal",
    'WBC' = "White-browed Coucal")

junk0 <- 
  data.frame(species = c("BC", "WBC"),
             stage = "HSR",
             difference = c(NA,NA))

sex_diff_survival_output2 <- 
  bind_rows(sex_diff_survival_output, junk0) %>% 
  mutate(difference = as.numeric(difference),
         stage = factor(stage, 
                        levels = c("HSR", "Nestling", "Groundling", 
                                   "Fledgling", "Adult")),
         species = factor(species, 
                          levels = c("BC", "WBC")))

# Figure 2a: plot the sex-biases in survival across the three stages
vital_rate_theme <- 
  theme_bw() +
  theme(
    text = element_text(family = "Verdana"),
    legend.position = "none",
    panel.background = element_rect(fill = "transparent", colour = NA),
    plot.background = element_rect(fill = "transparent", colour = NA),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.ticks.length = unit(0.1, "cm"),
    panel.border = element_blank(),
    panel.spacing = unit(0.3, "lines"),
    strip.background = element_blank()
  )

surv_diff_plot <-
  ggplot(aes(y = difference, x = stage, fill = stage), 
         data = slice_sample(group_by(sex_diff_survival_output, species, stage), n = 1000)) +
  geom_violin(draw_quantiles = c(0.025, 0.5, 0.975), color = "grey10",
              scale = "width", trim = TRUE, adjust = 1, size = 0.25) +
  vital_rate_theme +
  facet_grid(. ~ species, labeller = as_labeller(species_names)) +
  theme(axis.title.x = element_text(size = 7, colour = "black"),
        axis.text.x  = element_text(size = 6, angle = 45, hjust = 1, vjust = 1, 
                                    colour = "black"),
        axis.ticks.x = element_line(size = 0.2),
        axis.title.y = element_text(size = 7, hjust = 0.5, vjust = -3),
        axis.text.y  = element_text(size = 6, colour = "black"),
        axis.ticks.y = element_line(size = 0.2),
        plot.margin = unit(c(0.2, 1.35, 0.42, 0.3), "cm"),
        strip.text = element_text(size = 6, colour = "black", face = "italic")
  ) +
  scale_fill_manual(values = cbPalette) +
  scale_y_continuous(limits = c(-0.65, 0.65),
                     breaks = c(-0.5, -0.25, 0, 0.25, 0.5),
                     expand = c(0, 0), position = "left",
                     labels = c("-0.50","-0.25",
                                expression(phantom("-")*"0.00"),
                                expression(phantom("-")*"0.25"),
                                expression(phantom("-")*"0.50"))) +
  xlab("Life stage") +
  ylab("Sex bias\n") +
  scale_x_discrete(labels = c(
    "HSR" = "hatching\nsex ratio",
    "Nestling" = "nestling\nsurvival",
    "Groundling" = "groundling\nsurvival",
    "Fledgling" = "fledgling\nsurvival",
    "Adult" = "adult\nsurvival"))

background <-
  ggplot(aes(y = difference, x = stage, fill = stage),
         data = sex_diff_survival_output) +
  annotate("rect", xmin = -Inf, xmax = Inf, 
           ymin = -Inf, ymax = 0, alpha = 0.7,
           fill = pull(ggthemes_data$wsj$palettes$colors6[3,2])) +
  annotate("rect", xmin = -Inf, xmax = Inf, 
           ymin = 0, ymax = Inf, alpha = 0.7,
           fill = pull(ggthemes_data$wsj$palettes$colors6[2,2])) +
  annotate("text", x = c(5), y = c(-0.5),
           label = c("\u2640"), size = 4,
           family = "Menlo",
           vjust = c(0), hjust = c(0.5), colour = "grey10") +
  annotate("text", x = c(1), y = c(0.5),
           label = c("\u2642"), size = 4,
           family = "Menlo",
           vjust = c(1), hjust = c(0.5), colour = "grey10") +
  facet_grid(. ~ species, labeller = as_labeller(species_names)) +
  vital_rate_theme +
  theme(axis.title.x = element_text(size = 7, colour = "white"),
        axis.text.x  = element_text(size = 6, angle = 45, hjust = 1, 
                                    vjust = 1, colour = "white"),
        axis.ticks.x = element_line(size = 0.2, colour = "white"),
        axis.title.y = element_text(size = 7, hjust = 0.5, vjust = -1, 
                                    colour = "white"),
        axis.text.y  = element_text(size = 6, colour = "white"),
        axis.ticks.y = element_line(size = 0.2, colour = "white"),
        plot.margin = unit(c(0.2, 1.35, 1.55, 0.3), "cm"),
        strip.text = element_text(size = 6, colour = "white")
  ) +
  scale_x_discrete(labels = c(
    "HSR" = "hatching\nsex ratio",
    "Nestling" = "nestling\nsurvival",
    "Groundling" = "groundling\nsurvival",
    "Fledgling" = "fledgling\nsurvival",
    "Adult" = "adult\nsurvival")) +
  scale_y_continuous(limits = c(-0.65, 0.65), 
                     breaks = c(-0.5, -0.25, 0, 0.25, 0.5), 
                     expand = c(0, 0), position = "left",
                     labels = c("-0.50","-0.25", 
                                expression(phantom("-")*"0.00"), 
                                expression(phantom("-")*"0.25"),
                                expression(phantom("-")*"0.50"))) +
  xlab("Life stage") +
  ylab("Sex bias\n")

Derived Adult Sex Ratio

ASR_boot <- 
  bind_rows(BC_hazard_rate_boot_tidy$ASR_ests_boot,
            WBC_hazard_rate_boot_tidy$ASR_ests_boot) %>%
  mutate(species = factor(species, levels = c("BC", "WBC")))

CI <- 0.95

ASR_boot_summary <- 
  ASR_boot %>%
  dplyr::group_by(species) %>%
  dplyr::summarise(ucl_ASR = stats::quantile(M_Adult, (1 - CI)/2, na.rm = TRUE),
                   lcl_ASR = stats::quantile(M_Adult, 1 - (1 - CI)/2, na.rm = TRUE),
                   avg_ASR = mean(M_Adult),
                   med_ASR = median(M_Adult),
                   max_ASR = max(M_Adult),
                   min_ASR = min(M_Adult))

Figure_2b <- 
  ggplot() +
  annotate("rect", xmin=0, xmax=0.5, ymin=0, ymax=610, alpha=0.6,
           fill= pull(ggthemes_data$wsj$palettes$colors6[3,2])) +
  annotate("rect", xmin=0.5, xmax=1, ymin=0, ymax=610, alpha=0.6,
           fill= pull(ggthemes_data$wsj$palettes$colors6[2,2])) +
  annotate("text", x = c(0.05), y = c(305),
           label = c("\u2640"), size = 4, colour = "grey10",
           family="Menlo", vjust = c(0), hjust = c(0.5)) +
  annotate("text", x = c(0.95), y = c(305),
           label = c("\u2642"), size = 4, colour = "grey10",
           family="Menlo", vjust = c(1), hjust = c(0.5)) +
  geom_histogram(binwidth = 0.01, data = ASR_boot, aes(x = M_Adult), fill = "grey30") +
  geom_errorbarh(data = ASR_boot_summary, aes(y = 600, x = lcl_ASR, xmin = lcl_ASR, xmax = ucl_ASR), 
                 color = "black", size = 0.3, linetype = "solid") +
  coord_flip() +
  facet_grid(. ~ species) +
  theme_bw() +
  theme(text = element_text(family="Verdana"),
        legend.position = "none",
        axis.title.x = element_text(size=7, vjust=-0.1),
        axis.text.x  = element_text(size=6, angle = 45, hjust = 1, colour = "black"),
        axis.title.y = element_text(size=7, hjust=0.5, vjust = -2.75, margin = margin(0, 11, 0, 0)),
        axis.text.y  = element_text(size=6, colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.y = element_line(size = 0.2, colour = "grey40"),
        axis.ticks.length = unit(0.1, "cm"),
        axis.ticks.x = element_line(size = 0.2, colour = "grey40"),
        panel.border = element_blank(),
        plot.margin = unit(c(0.2, 1.35, 0.405, 0.08), "cm"),
        panel.spacing = unit(0.3, "lines"),
        strip.background = element_blank(),
        strip.text = element_blank()) +
  ylab("Bootstrap frequency") +
  xlab("Adult sex ratio\n(proportion \u2642, 95% CI)") +
  scale_x_continuous(limits = c(0, 1), expand = c(0, 0)) +
  scale_y_continuous(limits = c(0, 610), expand = c(0, 0), breaks=c(0, 100, 200, 300, 400, 500))

Life-table Response Experiment

Functions for the LTRE analysis

sensitivity_analysis()

Sensitivity analysis function takes the vital rate summary of the bootstrap procedure and conducts a perturbation analysis on each rate to assess how a proportional change in a given vital rate changes the ASR arguments: - vital_rate_summary: - matrix_str: - h: - k: - HSR: - niter: - ASR: - lambda:

sensitivity_analysis <-
  function(vital_rate_summary, matrix_str, h = 1, k = 4, 
           HSR, 
           niter = 1000, ASR, lambda){
    
    # make a list of all parameters
    vr <-
      list(F_Nestling_survival = vital_rate_summary$F_Nestling_survival,
           F_Groundling_survival = vital_rate_summary$F_Groundling_survival,
           F_Fledgling_survival = vital_rate_summary$F_Fledgling_survival,
           F_Adult_survival = vital_rate_summary$F_Adult_survival,
           M_Nestling_survival = vital_rate_summary$M_Nestling_survival,
           M_Groundling_survival = vital_rate_summary$M_Groundling_survival,
           M_Fledgling_survival = vital_rate_summary$M_Fledgling_survival,
           M_Adult_survival = vital_rate_summary$M_Adult_survival)
    
    # number of stages in the matrix
    no_stages <- sqrt(length(matrix_str))
    
    # Define plover life-stages of the Ceuta snowy plover matrix model
    stages <- c("F_1st_year",  "F_Adult",  "M_1st_year",  "M_Adult")
    
    # an empty t by x matrix
    stage <- matrix(numeric(no_stages * niter), nrow = no_stages)
    
    # an empty t vector to store the population sizes
    pop <- numeric(niter)
    
    # dataframe to store the perturbation results
    ASR_pert_results <-
      data.frame(parameter = c("F_Nestling_survival", "F_Groundling_survival", 
                               "F_Fledgling_survival", "F_Adult_survival",
                               "M_Nestling_survival", "F_Groundling_survival", 
                               "M_Fledgling_survival", "M_Adult_survival",
                               "h", "k", "HSR", "ISR"),
                 sensitivities = numeric(length(vr) + 4),
                 elasticities = numeric(length(vr) + 4))
    
    lambda_pert_results <-
      data.frame(parameter = c("F_Nestling_survival", "F_Groundling_survival", 
                               "F_Fledgling_survival", "F_Adult_survival",
                               "M_Nestling_survival", "F_Groundling_survival", 
                               "M_Fledgling_survival", "M_Adult_survival",
                               "h", "k", "HSR", "ISR"),
                 sensitivities = numeric(length(vr) + 4),
                 elasticities = numeric(length(vr) + 4))
    
    # specifiy how many survival rates there are
    n <- length(vr)
    
    # create vectors of perturbations to test on parameters of the matrix model
    vr_nums <- seq(0, 1, 0.01) # proportional changes in survival and HSR (i.e., between 0 an 1)
    h_nums <- seq(0, 2, 0.02) # proportional changes in h index (i.e., between 0 and 2)
    k_nums <- seq(3, 5, 0.02) # proportional changes in k (i.e, between 3 and 5)
    
    # create empty dataframes to store the perturbation results for ASR and lambda
    vr_pert_ASR <- matrix(numeric(n * length(vr_nums)),
                          ncol = n, dimnames = list(vr_nums, names(vr)))
    
    h_pert_ASR <- matrix(numeric(length(h_nums)),
                         ncol = 1, dimnames = list(h_nums, "h"))
    
    k_pert_ASR <- matrix(numeric(length(k_nums)),
                         ncol = 1, dimnames = list(k_nums, "k"))
    
    HSR_pert_ASR <- matrix(numeric(length(vr_nums)),
                           ncol = 1, dimnames = list(vr_nums, "HSR"))
    
    vr_pert_lambda <- matrix(numeric(n * length(vr_nums)),
                             ncol = n, dimnames = list(vr_nums, names(vr)))
    
    h_pert_lambda <- matrix(numeric(length(h_nums)),
                            ncol = 1, dimnames = list(h_nums, "h"))
    
    k_pert_lambda <- matrix(numeric(length(k_nums)),
                            ncol = 1, dimnames = list(k_nums, "k"))
    
    HSR_pert_lambda <- matrix(numeric(length(vr_nums)),
                              ncol = 1, dimnames = list(vr_nums, "HSR"))
    
    #### perturbation of survival rates ####
    for (g in 1:n) { # pick a column (i.e., a variable) 
      vr2 <- vr # reset the vital rates to the original
      
      for (i in 1:length(vr_nums)) # pick a perturbation level
      {
        
        vr2[[g]] <- vr_nums[i] # specify the vital rate with the new perturbation level
        
        A <- matrix(sapply(matrix_str, eval, vr2, NULL), 
                    nrow = sqrt(length(matrix_str)), byrow=TRUE, 
                    dimnames = list(stages, stages)) # build the matrix with the new value
        
        # reset the starting stage distribution for simulation (all with 10 individuals)
        m <- rep(10, no_stages) 
        
        for (j in 1:niter) { # project the matrix through t iteration
          # stage distribution at time t
          stage[,j] <- m
          
          # population size at time t
          pop[j] <- sum(m)
          
          # number of male adults at time t
          M2 <- stage[4, j] * A["M_Adult", "M_Adult"]
          
          # number of female adults at time t
          F2 <- stage[2, j] * A["F_Adult", "F_Adult"]
          
          # Female freq-dep fecundity of Female chicks
          A["F_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * (1 - HSR) )
          
          # Female freq-dep fecundity of Male chicks
          A["M_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * HSR)
          
          # Male freq-dep fecundity of Female chicks
          A["F_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * (1 - HSR))
          
          # Male freq-dep fecundity of Male chicks
          A["M_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * HSR)
          
          # define the new n (i.e., new stage distribution at time t)
          m <- A %*% (m)
        }
        
        # define rownames of stage matrix
        rownames(stage) <- rownames(A)
        
        # define colnames of stage matrix
        colnames(stage) <- 0:(niter - 1)
        
        # calculate the proportional stable stage distribution
        stage <- apply(stage, 2, function(x) x / sum(x))
        
        # define stable stage as the last stage
        stable.stage <- stage[, niter]
        
        # calc ASR as the proportion of the adult stable stage class that is male
        vr_pert_ASR[i, g] <- stable.stage[no_stages] / (stable.stage[no_stages/2] + 
                                                          stable.stage[no_stages])
        
        # calc lambda as the pop change in the counts of the last two iterations
        vr_pert_lambda[i, g] <- pop[niter]/pop[niter - 1]
      }
      
      # get the spline function of ASR
      spl_ASR <- 
        smooth.spline(vr_pert_ASR[which(!is.na(vr_pert_ASR[, g])), g] ~ 
                        names(vr_pert_ASR[which(!is.na(vr_pert_ASR[, g])), g]))
      
      # estimate the slope of the tangent of the spline at the vital rate
      ASR_pert_results[g, 2] <- predict(spl_ASR, x = vr[[g]], deriv = 1)$y
      
      # re-scale sensitivity into elasticity
      ASR_pert_results[g, 3] <- vr[[g]] / ASR * ASR_pert_results[g, 2]
      
      # do the same steps but for lambda
      spl_lambda <- 
        smooth.spline(vr_pert_lambda[which(!is.na(vr_pert_lambda[, g])), g] ~ 
                        names(vr_pert_lambda[which(!is.na(vr_pert_lambda[, g])), g]))
      
      lambda_pert_results[g, 2] <- predict(spl_lambda, x = vr[[g]], deriv = 1)$y
      
      lambda_pert_results[g, 3] <- vr[[g]] / lambda * lambda_pert_results[g, 2]
    }
    
    #### perturbation of the h index parameter ####
    for (i in 1:length(h_nums)) # pick a perturbation level
    {
      A <- matrix(sapply(matrix_str, eval, vr, NULL), 
                  nrow = sqrt(length(matrix_str)), byrow=TRUE, 
                  dimnames = list(stages, stages)) 
      
      # reset the starting stage distribution for simulation (all with 10 individuals)
      m <- rep(10, no_stages) 
      
      for (j in 1:niter) { # project the matrix through t iteration
        # stage distribution at time t
        stage[,j] <- m
        
        # population size at time t
        pop[j] <- sum(m)
        
        # number of male adults at time t
        M2 <- stage[4, j] * A["M_Adult", "M_Adult"]
        
        # number of female adults at time t
        F2 <- stage[2, j] * A["F_Adult", "F_Adult"]
        
        # Female freq-dep fecundity of Female chicks
        A["F_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h_nums[i])) * (1 - HSR) )
        
        # Female freq-dep fecundity of Male chicks
        A["M_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h_nums[i])) * HSR)
        
        # Male freq-dep fecundity of Female chicks
        A["F_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h_nums[i])) * (1 - HSR))
        
        # Male freq-dep fecundity of Male chicks
        A["M_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h_nums[i])) * HSR)
        
        # define the new n (i.e., new stage distribution at time t)
        m <- A %*% (m)
      }
      # define rownames of stage matrix
      rownames(stage) <- rownames(A)
      
      # define colnames of stage matrix
      colnames(stage) <- 0:(niter - 1)
      
      # calculate the proportional stable stage distribution
      stage <- apply(stage, 2, function(x) x / sum(x))
      
      # define stable stage as the last stage
      stable.stage <- stage[, niter]
      
      # calc ASR as the proportion of the adult stable stage class that is male
      h_pert_ASR[i,] <- stable.stage[no_stages] / (stable.stage[no_stages / 2] + 
                                                     stable.stage[no_stages])
      
      # calc lambda as the pop change in the counts of the last two iterations
      h_pert_lambda[i, ] <- pop[niter]/pop[niter - 1]
      
    }
    # get the spline function of ASR
    spl_ASR <- 
      smooth.spline(h_pert_ASR[which(!is.na(h_pert_ASR)), 1] ~ 
                      names(h_pert_ASR[which(!is.na(h_pert_ASR)), ]))
    
    # estimate the slope of the tangent of the spline at the vital rate
    ASR_pert_results[n + 1, 2] <- predict(spl_ASR, x = h, deriv = 1)$y
    
    # re-scale sensitivity into elasticity
    ASR_pert_results[n + 1, 3] <- h / ASR * ASR_pert_results[n + 1, 2]
    
    # do the same steps but for lambda
    spl_lambda <- 
      smooth.spline(h_pert_lambda[which(!is.na(h_pert_lambda)), 1] ~ 
                      names(h_pert_lambda[which(!is.na(h_pert_lambda)), ]))
    lambda_pert_results[n + 1, 2] <- predict(spl_lambda, x = h, deriv = 1)$y
    lambda_pert_results[n + 1, 3] <- h / lambda * lambda_pert_results[n + 1, 2]
    
    #### perturbation of k parameter ####
    for (i in 1:length(k_nums)) # pick a perturbation level
    {
      A <- matrix(sapply(matrix_str, eval, vr, NULL), 
                  nrow = sqrt(length(matrix_str)), byrow=TRUE, 
                  dimnames = list(stages, stages))
      
      # reset the starting stage distribution for simulation (all with 10 individuals)
      m <- rep(10, no_stages) 
      
      for (j in 1:niter) { # project the matrix through t iteration
        # stage distribution at time t
        stage[,j] <- m
        
        # population size at time t
        pop[j] <- sum(m)
        
        # number of male adults at time t
        M2 <- stage[4, j] * A["M_Adult", "M_Adult"]
        
        # number of female adults at time t
        F2 <- stage[2, j] * A["F_Adult", "F_Adult"]
        
        # Female freq-dep fecundity of Female chicks
        A["F_1st_year", "F_Adult"] <- ((k_nums[i] * M2) / (M2 + (F2 / h)) * (1 - HSR) )
        
        # Female freq-dep fecundity of Male chicks
        A["M_1st_year", "F_Adult"] <- ((k_nums[i] * M2) / (M2 + (F2 / h)) * HSR)
        
        # Male freq-dep fecundity of Female chicks
        A["F_1st_year", "M_Adult"] <- ((k_nums[i] * F2) / (M2 + (F2 / h)) * (1 - HSR))
        
        # Male freq-dep fecundity of Male chicks
        A["M_1st_year", "M_Adult"] <- ((k_nums[i] * F2) / (M2 + (F2 / h)) * HSR)
        
        # define the new n (i.e., new stage distribution at time t)
        m <- A %*% (m)
      }
      # define rownames of stage matrix
      rownames(stage) <- rownames(A)
      
      # define colnames of stage matrix
      colnames(stage) <- 0:(niter - 1)
      
      # calculate the proportional stable stage distribution
      stage <- apply(stage, 2, function(x) x / sum(x))
      
      # define stable stage as the last stage
      stable.stage <- stage[, niter]
      
      # calc ASR as the proportion of the adult stable stage class that is male
      k_pert_ASR[i,] <- stable.stage[no_stages] / (stable.stage[no_stages/2] + 
                                                     stable.stage[no_stages])
      
      # calc lambda as the pop change in the counts of the last two iterations
      k_pert_lambda[i, ] <- pop[niter]/pop[niter - 1]
    }
    
    # get the spline function of ASR
    spl_ASR <- 
      smooth.spline(k_pert_ASR[which(!is.na(k_pert_ASR)), 1] ~ 
                      names(k_pert_ASR[which(!is.na(k_pert_ASR)), ]))
    
    # estimate the slope of the tangent of the spline at the vital rate
    
    ASR_pert_results[n + 2, 2] <- predict(spl_ASR, x = k, deriv = 1)$y
    
    # re-scale sensitivity into elasticity
    ASR_pert_results[n + 2, 3] <- k / ASR * ASR_pert_results[n + 2, 2]
    
    # do the same steps but for lambda
    spl_lambda <- 
      smooth.spline(h_pert_lambda[which(!is.na(h_pert_lambda)), 1] ~ 
                      names(h_pert_lambda[which(!is.na(h_pert_lambda)), ]))
    lambda_pert_results[n + 2, 2] <- predict(spl_lambda, x = k, deriv = 1)$y
    lambda_pert_results[n + 2, 3] <- k / lambda * lambda_pert_results[n + 2, 2]
    
    #### perturbation of HSR ####
    for (i in 1:length(vr_nums)) # pick a perturbation level
    {
      A <- matrix(sapply(matrix_str, eval, vr, NULL), 
                  nrow = sqrt(length(matrix_str)), byrow=TRUE, 
                  dimnames = list(stages, stages))
      
      # reset the starting stage distribution for simulation (all with 10 individuals)
      m <- rep(10, no_stages) 
      
      for (j in 1:niter) { # project the matrix through t iteration
        # stage distribution at time t
        stage[,j] <- m
        
        # population size at time t
        pop[j] <- sum(m)
        
        # number of male adults at time t
        M2 <- stage[4, j] * A["M_Adult", "M_Adult"]
        
        # number of female adults at time t
        F2 <- stage[2, j] * A["F_Adult", "F_Adult"]
        
        # Female freq-dep fecundity of Female chicks
        A["F_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * (1 - vr_nums[i]) )
        
        # Female freq-dep fecundity of Male chicks
        A["M_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * vr_nums[i])
        
        # Male freq-dep fecundity of Female chicks
        A["F_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * (1 - vr_nums[i]))
        
        # Male freq-dep fecundity of Male chicks
        A["M_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * vr_nums[i])
        
        # define the new n (i.e., new stage distribution at time t)
        m <- A %*% (m)
        
      }
      # define rownames of stage matrix
      rownames(stage) <- rownames(A)
      
      # define colnames of stage matrix
      colnames(stage) <- 0:(niter - 1)
      
      # calculate the proportional stable stage distribution
      stage <- apply(stage, 2, function(x) x / sum(x))
      
      # define stable stage as the last stage
      stable.stage <- stage[, niter]
      
      # calc ASR as the proportion of the adult stable stage class that is male
      HSR_pert_ASR[i,] <- stable.stage[no_stages] / (stable.stage[no_stages/2] + 
                                                       stable.stage[no_stages])
      
      # calc lambda as the pop change in the counts of the last two iterations
      HSR_pert_lambda[i, ] <- pop[niter] / pop[niter - 1]
      
    }
    # get the spline function of ASR
    spl_ASR <- 
      smooth.spline(HSR_pert_ASR[which(!is.na(HSR_pert_ASR)), 1] ~ 
                      names(HSR_pert_ASR[which(!is.na(HSR_pert_ASR)), ]))
    
    # estimate the slope of the tangent of the spline at the vital rate    
    ASR_pert_results[n + 3, 2] <- predict(spl_ASR, x = HSR, deriv = 1)$y
    
    # re-scale sensitivity into elasticity
    ASR_pert_results[n + 3, 3] <- HSR / ASR * ASR_pert_results[n + 3, 2]
    
    # do the same steps but for lambda
    spl_lambda <- 
      smooth.spline(h_pert_lambda[which(!is.na(h_pert_lambda)), 1] ~ 
                      names(h_pert_lambda[which(!is.na(h_pert_lambda)), ]))
    lambda_pert_results[n + 3, 2] <- predict(spl_lambda, x = HSR, deriv = 1)$y
    lambda_pert_results[n + 3, 3] <- HSR/lambda * lambda_pert_results[n + 3, 2]
    
    #### store all results into a list ----
    result <- list(ASR_pert_results = ASR_pert_results,
                   lambda_pert_results = lambda_pert_results,
                   HSR_pert_ASR = HSR_pert_ASR,
                   HSR_pert_lambda = HSR_pert_lambda,
                   k_pert_ASR = k_pert_ASR,
                   k_pert_lambda = k_pert_lambda,
                   h_pert_ASR = h_pert_ASR,
                   h_pert_lambda = h_pert_lambda,
                   vr_pert_ASR = vr_pert_ASR,
                   vr_pert_lambda = vr_pert_lambda)
    
  }

matrix_structure

matrix_structure is an expression that determines the structure of the two-sex matrix from a set of vital rates

matrix_structure <- expression(
  # top row of matrix
  0, NA, 0, NA,
  
  # second row of matrix
  (F_Nestling_survival * F_Groundling_survival * F_Fledgling_survival), F_Adult_survival, 0, 0,
  
  # third row of matrix
  0, NA, 0, NA,
  
  # fourth row of matrix
  0, 0, (M_Nestling_survival * M_Groundling_survival * M_Fledgling_survival), M_Adult_survival
  
)

make_treat_matrix()

make_treat_matrix() fucntion makes a treatment matrix from a summary of the bootstrapped vital rates

arguments:

  • survival_rates_boot_summary:
  • species_name:
  • h:
  • k:
  • HSR:
make_treat_matrix <- 
  function(survival_rates_boot_summary, 
           species_name, 
           h, 
           k, 
           HSR){
    
    list(F_Nestling_survival = filter(survival_rates_boot_summary,
                                      species == species_name & vital_rate == "Female_nestling_survival")[, 4],
         F_Groundling_survival = filter(survival_rates_boot_summary,
                                        species == species_name & vital_rate == "Female_groundling_survival")[, 4],
         F_Fledgling_survival = filter(survival_rates_boot_summary,
                                       species == species_name & vital_rate == "Female_fledgling_survival")[, 4],
         F_Adult_survival = filter(survival_rates_boot_summary,
                                   species == species_name & vital_rate == "Female_adult_survival")[, 4],
         M_Nestling_survival = filter(survival_rates_boot_summary,
                                      species == species_name & vital_rate == "Male_nestling_survival")[, 4],
         M_Groundling_survival = filter(survival_rates_boot_summary,
                                        species == species_name & vital_rate == "Male_groundling_survival")[, 4],
         M_Fledgling_survival = filter(survival_rates_boot_summary,
                                       species == species_name & vital_rate == "Male_fledgling_survival")[, 4],
         M_Adult_survival = filter(survival_rates_boot_summary,
                                   species == species_name & vital_rate == "Male_adult_survival")[, 4],
         Egg_survival = filter(survival_rates_boot_summary,
                               species == species_name & vital_rate == "NA_egg_survival")[, 4],
         
         # Define h (harem size, h = 1 is monogamy) and k (clutch size)
         h = h,
         k = k,
         
         # Define primary sex ratio
         HSR = HSR)
  }

make_mprime_matrix()

make_mprime_matrix() fucntion makes a prime-matrix (i.e., a matrix halfway between the treatment matrix and the unbiased matrix) from a set of vital rate summaries

arguments:

  • survival_rates_boot_summary:
  • species_name:
  • h:
  • k:
  • HSR:
  • sex:
make_mprime_matrix <- 
  function(survival_rates_boot_summary, 
           species_name, 
           h, 
           k, 
           HSR, 
           sex){
    if(sex == "male"){
      list(F_Nestling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_nestling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_nestling_survival")[, 4]) / 2,
           F_Groundling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_groundling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_groundling_survival")[, 4]) / 2,
           F_Fledgling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_fledgling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_fledgling_survival")[, 4]) / 2,
           F_Adult_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_adult_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_adult_survival")[, 4]) / 2,
           M_Nestling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Male_nestling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_nestling_survival")[, 4]) / 2,
           M_Groundling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Male_groundling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_groundling_survival")[, 4]) / 2,
           M_Fledgling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Male_fledgling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_fledgling_survival")[, 4]) / 2,
           M_Adult_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Male_adult_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_adult_survival")[, 4]) / 2,
           Egg_survival = filter(survival_rates_boot_summary,
                                 species == species_name)[15, 4],
           
           # Define h (harem size, h = 1 is monogamy) and k (clutch size)
           h = (h + 1) / 2,
           k = k,
      
           # Define primary sex ratio
           HSR = (HSR + 0.5) / 2)
           

    }
    else{
      list(F_Nestling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_nestling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Female_nestling_survival")[, 4]) / 2,
           F_Groundling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_groundling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Female_groundling_survival")[, 4]) / 2,
           F_Fledgling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_fledgling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Female_fledgling_survival")[, 4]) / 2,
           F_Adult_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_adult_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Female_adult_survival")[, 4]) / 2,
           M_Nestling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_nestling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_nestling_survival")[, 4]) / 2,
           M_Groundling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_groundling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_groundling_survival")[, 4]) / 2,
           M_Fledgling_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_fledgling_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_fledgling_survival")[, 4]) / 2,
           M_Adult_survival = (filter(survival_rates_boot_summary,
                                         species == species_name & vital_rate == "Female_adult_survival")[, 4] +
                                    filter(survival_rates_boot_summary,
                                           species == species_name & vital_rate == "Male_adult_survival")[, 4]) / 2,
           Egg_survival = filter(survival_rates_boot_summary,
                                 species == species_name)[15, 4],
           
           # Define h (harem size, h = 1 is monogamy) and k (clutch size)
           h = (h + 1) / 2,
           k = k,
           
           # Define primary sex ratio
           HSR = (HSR + 0.5) / 2)
    }
  }

LTRE_analysis()

LTRE_analysis() compares the relative difference in ASR response of an unbiased two-sex matrix vs. the observed sex-specific matrix

arguments:

  • Mprime_sens:
  • matrix_str:
  • vital_rates:
  • species_name:
  • sex:
LTRE_analysis <-
  function(Mprime_sens, matrix_str, vital_rates, species_name, sex){
    
    # make empty dataframes to store LTRE results for ASR and lambda
    LTRE_ASR <-
      data.frame(parameter = c("Nestling survival", "Groundling survival",
                               "Fledgling survival", "Adult survival",
                               "Hatching sex ratio", 
                               "Mating system"),
                 contribution = numeric(6))
    LTRE_lambda <-
      data.frame(parameter = c("Nestling survival", "Groundling survival",
                               "Fledgling survival", "Adult survival",
                               "Hatching sex ratio",
                               "Mating system"),
                 contribution = numeric(6))
    
    # run a for loop to extract the parameter contributions
    if (sex == "male"){
      # male rates scenario
      for(i in 1:nrow(LTRE_ASR))
      {
        LTRE_ASR[i, 2] <-
          
          # survival rates
          ifelse(i < 5, (vital_rates[[i + 4]] - vital_rates[[i]]) * 
                   Mprime_sens$ASR_pert_results$sensitivities[i + 4],
                 
                 # HSR
                 ifelse(i == 5, (vital_rates[[12]] - 0.5) * 
                          Mprime_sens$ASR_pert_results$sensitivities[11],
                               
                               # mating system
                               (1 - vital_rates[[10]]) * 
                                 Mprime_sens$ASR_pert_results$sensitivities[9]))
      }
      for(i in 1:nrow(LTRE_lambda))
      {
        # survival rates
        ifelse(i < 5, (vital_rates[[i + 4]] - vital_rates[[i]]) * 
                 Mprime_sens$lambda_pert_results$sensitivities[i + 4],
               
               # HSR
               ifelse(i == 5, (vital_rates[[12]] - 0.5) * 
                        Mprime_sens$lambda_pert_results$sensitivities[11],
                             
                             # mating system
                             (1 - vital_rates[[10]]) * 
                               Mprime_sens$lambda_pert_results$sensitivities[9]))
      }
    }
    else{
      # female rates scenario
      for(i in 1:nrow(LTRE_ASR))
      {
        LTRE_ASR[i, 2] <-
          
          # survival rates
          ifelse(i < 5, (vital_rates[[i]] - vital_rates[[i + 4]]) * 
                   Mprime_sens$ASR_pert_results$sensitivities[i],
                 
                 # HSR
                 ifelse(i == 5, (vital_rates[[12]] - 0.5) * 
                          Mprime_sens$ASR_pert_results$sensitivities[11],
                        
                               # mating system
                               (1 - vital_rates[[10]]) * 
                                 Mprime_sens$ASR_pert_results$sensitivities[9]))
      }
      
      for(i in 1:nrow(LTRE_lambda))
      {
        # survival rates
        ifelse(i < 5, (vital_rates[[i]] - vital_rates[[i + 4]]) * 
                 Mprime_sens$lambda_pert_results$sensitivities[i],
               
               # HSR
               ifelse(i == 5, (vital_rates[[12]] - 0.5) * 
                        Mprime_sens$lambda_pert_results$sensitivities[11],
                             
                             # mating system
                             (1 - vital_rates[[10]]) * 
                               Mprime_sens$lambda_pert_results$sensitivities[9]))
      }
    }
    LTRE_ASR$parameter <- 
      factor(LTRE_ASR$parameter, levels = c("Adult survival",
                                            "Fledgling survival",
                                            "Groundling survival",
                                            "Nestling survival",
                                            "Hatching sex ratio",
                                            "Mating system"))
    LTRE_lambda$parameter <- 
      factor(LTRE_lambda$parameter, levels = c("Adult survival",
                                               "Fledgling survival",
                                               "Groundling survival",
                                               "Nestling survival",
                                               "Hatching sex ratio",
                                               "Mating system"))
    
    LTRE_ASR$model <- "ASR"
    LTRE_ASR$species <- species_name
    
    LTRE_lambda$model <- "lambda"
    LTRE_lambda$species <- species_name
    
    LTRE_results <- list(LTRE_ASR = LTRE_ASR,
                         LTRE_lambda = LTRE_lambda)
    
    LTRE_results
  }

LTRE_contributions_check()

LTRE_contributions_check() takes the results of the LTRE analysis and checks the total contribution that each vital rate has on ASR. In general, the sum of the LTRE contribution should roughly equal the difference in ASR between the observed scenario and the unbiased scenario.

arguments:

  • M_matrix_vital_rates:
  • Mprime_sensitivities:
  • M_matrix_ASR:
  • scenario:
LTRE_contributions_check <- 
  function(M_matrix_vital_rates, Mprime_sensitivities, M_matrix_ASR, scenario){
    if (scenario == "male")
      contribution_sum <- 
        sum(
          (M_matrix_vital_rates[[1]] - M_matrix_vital_rates[[5]]) * 
            Mprime_sensitivities[1],
          
          (M_matrix_vital_rates[[2]] - M_matrix_vital_rates[[6]]) * 
            Mprime_sensitivities[2],
          
          (M_matrix_vital_rates[[3]] - M_matrix_vital_rates[[7]]) * 
            Mprime_sensitivities[3],
          
          (M_matrix_vital_rates[[4]] - M_matrix_vital_rates[[8]]) * 
            Mprime_sensitivities[4],
          
          (M_matrix_vital_rates[[5]] - M_matrix_vital_rates[[5]]) * 
            Mprime_sensitivities[1],
          
          (M_matrix_vital_rates[[6]] - M_matrix_vital_rates[[6]]) * 
            Mprime_sensitivities[2],
          
          (M_matrix_vital_rates[[7]] - M_matrix_vital_rates[[7]]) * 
            Mprime_sensitivities[3],
          
          (M_matrix_vital_rates[[8]] - M_matrix_vital_rates[[8]]) * 
            Mprime_sensitivities[4],
          
          (1 - M_matrix_vital_rates[[10]]) *
            Mprime_sensitivities[9],
          
          (M_matrix_vital_rates[[12]] - 0.5) * Mprime_sensitivities[11],
          
          (M_matrix_vital_rates[[13]] - 0.5) * Mprime_sensitivities[12]
          
        )
    
    else
      contribution_sum <- 
        sum(
          (M_matrix_vital_rates[[5]] - M_matrix_vital_rates[[1]]) * 
            Mprime_sensitivities[5],
          
          (M_matrix_vital_rates[[6]] - M_matrix_vital_rates[[2]]) * 
            Mprime_sensitivities[6],
          
          (M_matrix_vital_rates[[7]] - M_matrix_vital_rates[[3]]) * 
            Mprime_sensitivities[7],
          
          (M_matrix_vital_rates[[8]] - M_matrix_vital_rates[[4]]) * 
            Mprime_sensitivities[8],
          
          (M_matrix_vital_rates[[1]] - M_matrix_vital_rates[[1]]) * 
            Mprime_sensitivities[5],
          
          (M_matrix_vital_rates[[2]] - M_matrix_vital_rates[[2]]) * 
            Mprime_sensitivities[6],
          
          (M_matrix_vital_rates[[3]] - M_matrix_vital_rates[[3]]) * 
            Mprime_sensitivities[7],
          
          (M_matrix_vital_rates[[4]] - M_matrix_vital_rates[[4]]) * 
            Mprime_sensitivities[8],
          
          (1 - M_matrix_vital_rates[[10]]) *
            Mprime_sensitivities[9],
          
          (M_matrix_vital_rates[[12]] - 0.5) * Mprime_sensitivities[11],
          
          (M_matrix_vital_rates[[13]] - 0.5) * Mprime_sensitivities[12]
        )
    ASR_bias <- abs(M_matrix_ASR - 0.5)
    absolute_difference <- abs(ASR_bias) - abs(contribution_sum)
    
    return(list(contribution_sum = as.vector(contribution_sum), 
                ASR_bias = as.vector(ASR_bias), 
                absolute_difference = as.vector(absolute_difference)))
  }

LTRE: Analysis

BC_hazard_rate_boot_tidy$vital_rate_ests_boot$iter <- 
  as.factor(BC_hazard_rate_boot_tidy$vital_rate_ests_boot$iter)
WBC_hazard_rate_boot_tidy$vital_rate_ests_boot$iter <- 
  as.factor(WBC_hazard_rate_boot_tidy$vital_rate_ests_boot$iter)

# summarize the vital rates
survival_rates_boot_summary <-
  bind_rows(BC_hazard_rate_boot_tidy$vital_rate_ests_boot, WBC_hazard_rate_boot_tidy$vital_rate_ests_boot) %>% 
  mutate(vital_rate = paste(sex, stage, rate, sep = "_")) %>% 
  Rmisc::summarySE(.,
                   measurevar = "value",
                   groupvars = c("vital_rate", "species"),
                   conf.interval = 0.95) %>% 
  arrange(species, vital_rate)


BC_VR_treat <- 
  make_treat_matrix(survival_rates_boot_summary = survival_rates_boot_summary,
                    h = 1/pull(filter(parameter_distributions, species == "BC" & trait == "mating_system"), mean),
                    HSR = pull(filter(parameter_distributions, species == "BC" & trait == "hatching_sex_ratio"), mean),
                    k = pull(filter(parameter_distributions, species == "BC" & trait == "clutch_size"), mean),
                    species_name = "BC")

BC_VR_mprime_male <- 
  make_mprime_matrix(survival_rates_boot_summary = survival_rates_boot_summary,
                     h = 1/pull(filter(parameter_distributions, species == "BC" & trait == "mating_system"), mean),
                     HSR = pull(filter(parameter_distributions, species == "BC" & trait == "hatching_sex_ratio"), mean),
                     k = pull(filter(parameter_distributions, species == "BC" & trait == "clutch_size"), mean),
                     species_name = "BC",
                     sex = "male")

BC_VR_mprime_female <- 
  make_mprime_matrix(survival_rates_boot_summary = survival_rates_boot_summary,
                     h = 1/pull(filter(parameter_distributions, species == "BC" & trait == "mating_system"), mean),
                     HSR = pull(filter(parameter_distributions, species == "BC" & trait == "hatching_sex_ratio"), mean),
                     k = pull(filter(parameter_distributions, species == "BC" & trait == "clutch_size"), mean),
                     species_name = "BC",
                     sex = "female")
WBC_VR_treat <- 
  make_treat_matrix(survival_rates_boot_summary = survival_rates_boot_summary,
                    h = 1/pull(filter(parameter_distributions, species == "WBC" & trait == "mating_system"), mean),
                    HSR = pull(filter(parameter_distributions, species == "WBC" & trait == "hatching_sex_ratio"), mean),
                    k = pull(filter(parameter_distributions, species == "WBC" & trait == "clutch_size"), mean),
                    species_name = "WBC")

WBC_VR_mprime_male <- 
  make_mprime_matrix(survival_rates_boot_summary = survival_rates_boot_summary,
                     h = 1/pull(filter(parameter_distributions, species == "WBC" & trait == "mating_system"), mean),
                     HSR = pull(filter(parameter_distributions, species == "WBC" & trait == "hatching_sex_ratio"), mean),
                     k = pull(filter(parameter_distributions, species == "WBC" & trait == "clutch_size"), mean),
                     species_name = "WBC",
                     sex = "male")

WBC_VR_mprime_female <- 
  make_mprime_matrix(survival_rates_boot_summary = survival_rates_boot_summary,
                     h = 1/pull(filter(parameter_distributions, species == "WBC" & trait == "mating_system"), mean),
                     HSR = pull(filter(parameter_distributions, species == "WBC" & trait == "hatching_sex_ratio"), mean),
                     k = pull(filter(parameter_distributions, species == "WBC" & trait == "clutch_size"), mean),
                     species_name = "WBC",
                     sex = "female")

BC_treatment_matrix <- coucal_matrix(BC_VR_treat)
BC_M_prime_matrix_male <- coucal_matrix(BC_VR_mprime_male)
BC_M_prime_matrix_female <- coucal_matrix(BC_VR_mprime_female)

WBC_treatment_matrix <- coucal_matrix(WBC_VR_treat)
WBC_M_prime_matrix_male <- coucal_matrix(WBC_VR_mprime_male)
WBC_M_prime_matrix_female <- coucal_matrix(WBC_VR_mprime_female)

BC_treatment_ASR_analysis <- 
  matrix_ASR(M = BC_treatment_matrix,  
             h = BC_VR_treat$h, 
             k = BC_VR_treat$k,
             HSR = BC_VR_treat$HSR,
             iterations = 100, 
             num_boot = 1, 
             iter_add = 1,
             species = "BC")

BC_ASR_treat <- BC_treatment_ASR_analysis$ASR
BC_ASR_treat
  M_Adult 
0.6901214 
WBC_treatment_ASR_analysis <- 
  matrix_ASR(M = WBC_treatment_matrix,  
             h = WBC_VR_treat$h, 
             k = WBC_VR_treat$k,
             HSR = WBC_VR_treat$HSR,
             iterations = 100, 
             num_boot = 1, 
             iter_add = 1,
             species = "BC")

WBC_ASR_treat <- WBC_treatment_ASR_analysis$ASR
WBC_ASR_treat
  M_Adult 
0.4936118 
BC_M_prime_ASR_analysis_male <- 
  matrix_ASR(M = BC_M_prime_matrix_male, 
             h = BC_VR_mprime_male$h, 
             k = BC_VR_mprime_male$k,
             HSR = BC_VR_mprime_male$HSR,
             iterations = 100,
             num_boot = 1,
             iter_add = 1,
             species = "BC")

BC_ASR_mprime_male <- BC_M_prime_ASR_analysis_male$ASR
BC_ASR_mprime_male
  M_Adult 
0.5855401 
BC_M_prime_ASR_analysis_female <- 
  matrix_ASR(M = BC_M_prime_matrix_female, 
             h = BC_VR_mprime_female$h, 
             k = BC_VR_mprime_female$k,
             HSR = BC_VR_mprime_female$HSR,
             iterations = 100,
             num_boot = 1,
             iter_add = 1,
             species = "BC")

BC_ASR_mprime_female <- BC_M_prime_ASR_analysis_female$ASR
BC_ASR_mprime_female
  M_Adult 
0.5797847 
WBC_M_prime_ASR_analysis_male <- 
  matrix_ASR(M = WBC_M_prime_matrix_male, 
             h = WBC_VR_mprime_male$h, 
             k = WBC_VR_mprime_male$k,
             HSR = WBC_VR_mprime_male$HSR,
             iterations = 100,
             num_boot = 1,
             iter_add = 1,
             species = "BC")

WBC_ASR_mprime_male <- WBC_M_prime_ASR_analysis_male$ASR
WBC_ASR_mprime_male
  M_Adult 
0.4951691 
WBC_M_prime_ASR_analysis_female <- 
  matrix_ASR(M = WBC_M_prime_matrix_female, 
             h = WBC_VR_mprime_female$h, 
             k = WBC_VR_mprime_female$k,
             HSR = WBC_VR_mprime_female$HSR,
             iterations = 100,
             num_boot = 1,
             iter_add = 1,
             species = "BC")

WBC_ASR_mprime_female <- WBC_M_prime_ASR_analysis_female$ASR
WBC_ASR_mprime_female
  M_Adult 
0.4983036 
BC_lambda_treat <- BC_treatment_ASR_analysis$lambda
BC_lambda_treat
[1] 0.8686866
BC_lambda_mprime_male <- BC_M_prime_ASR_analysis_male$lambda
BC_lambda_mprime_male
[1] 1.008295
BC_lambda_mprime_female <- BC_M_prime_ASR_analysis_female$lambda
BC_lambda_mprime_female
[1] 0.8823659
WBC_lambda_treat <- WBC_treatment_ASR_analysis$lambda
WBC_lambda_treat
[1] 0.9134844
WBC_lambda_mprime_male <- WBC_M_prime_ASR_analysis_male$lambda
WBC_lambda_mprime_male
[1] 0.920894
WBC_lambda_mprime_female <- WBC_M_prime_ASR_analysis_female$lambda
WBC_lambda_mprime_female
[1] 0.9242015
BC_treat_sensitivity_analysis <- 
  sensitivity_analysis(vital_rate_summary = BC_VR_treat, 
                       matrix_str = matrix_structure, 
                       h = BC_VR_treat$h, 
                       k = BC_VR_treat$k, 
                       HSR = BC_VR_treat$HSR, 
                       niter = 100, 
                       ASR = BC_ASR_treat,
                       lambda = BC_lambda_treat)

BC_Mprime_sensitivity_analysis_male <- 
  sensitivity_analysis(vital_rate_summary = BC_VR_mprime_male, 
                       matrix_str = matrix_structure, 
                       h = BC_VR_mprime_male$h, 
                       k = BC_VR_mprime_male$k, 
                       HSR = BC_VR_mprime_male$HSR, 
                       niter = 100, 
                       ASR = BC_ASR_mprime_male,
                       lambda = BC_lambda_mprime_male)

BC_Mprime_sensitivity_analysis_female <- 
  sensitivity_analysis(vital_rate_summary = BC_VR_mprime_female, 
                       matrix_str = matrix_structure, 
                       h = BC_VR_mprime_female$h, 
                       k = BC_VR_mprime_female$k, 
                       HSR = BC_VR_mprime_female$HSR, 
                       niter = 100, 
                       ASR = BC_ASR_mprime_female,
                       lambda = BC_lambda_mprime_female) 

WBC_treat_sensitivity_analysis <- 
  sensitivity_analysis(vital_rate_summary = WBC_VR_treat, 
                       matrix_str = matrix_structure, 
                       h = WBC_VR_treat$h, 
                       k = WBC_VR_treat$k, 
                       HSR = WBC_VR_treat$HSR, 
                       niter = 100, 
                       ASR = WBC_ASR_treat,
                       lambda = WBC_lambda_treat)

WBC_Mprime_sensitivity_analysis_male <- 
  sensitivity_analysis(vital_rate_summary = WBC_VR_mprime_male, 
                       matrix_str = matrix_structure, 
                       h = WBC_VR_mprime_male$h, 
                       k = WBC_VR_mprime_male$k, 
                       HSR = WBC_VR_mprime_male$HSR, 
                       niter = 100, 
                       ASR = WBC_ASR_mprime_male,
                       lambda = WBC_lambda_mprime_male) 

WBC_Mprime_sensitivity_analysis_female <- 
  sensitivity_analysis(vital_rate_summary = WBC_VR_mprime_female, 
                       matrix_str = matrix_structure, 
                       h = WBC_VR_mprime_female$h, 
                       k = WBC_VR_mprime_female$k, 
                       HSR = WBC_VR_mprime_female$HSR, 
                       niter = 100, 
                       ASR = WBC_ASR_mprime_female,
                       lambda = WBC_lambda_mprime_female)
BC_LTRE_male <- 
  LTRE_analysis(Mprime_sens = BC_Mprime_sensitivity_analysis_male, 
                matrix_str = matrix_str, 
                vital_rates = BC_VR_treat,
                species_name = "BC",
                sex = "male")

BC_LTRE_female <- 
  LTRE_analysis(Mprime_sens = BC_Mprime_sensitivity_analysis_female, 
                matrix_str = matrix_str, 
                vital_rates = BC_VR_treat,
                species_name = "BC",
                sex = "female")

WBC_LTRE_male <- 
  LTRE_analysis(Mprime_sens = WBC_Mprime_sensitivity_analysis_male, 
                matrix_str = matrix_str, 
                vital_rates = WBC_VR_treat,
                species_name = "WBC",
                sex = "male")

WBC_LTRE_female <- 
  LTRE_analysis(Mprime_sens = WBC_Mprime_sensitivity_analysis_female, 
                matrix_str = matrix_str, 
                vital_rates = WBC_VR_treat,
                species_name = "WBC",
                sex = "female")

LTRE_coucal_male_ASR <- 
  bind_rows(BC_LTRE_male$LTRE_ASR, WBC_LTRE_male$LTRE_ASR) %>% 
  mutate(sex = "male")

LTRE_coucal_female_ASR <- 
  bind_rows(BC_LTRE_female$LTRE_ASR, WBC_LTRE_female$LTRE_ASR) %>% 
  mutate(sex = "female")

LTRE_coucal_ASR <- rbind(LTRE_coucal_male_ASR, LTRE_coucal_female_ASR)

LTRE_coucal_ASR$parameter <- 
  factor(LTRE_coucal_ASR$parameter, 
         levels = c("Hatching sex ratio",
                    "Nestling survival",
                    "Groundling survival",
                    "Fledgling survival",
                    "Adult survival",
                    "Mating system"))

LTRE Results

# custom color palette for the plotting of Juvenile and Adult stats
cbPalette <- c("#D9D9D9", "#D9D9D9", "#D9D9D9", 
               "#D9D9D9", "#A6A6A6", "#A6A6A6",
               "#A6A6A6")

species_names <- 
  c('BC' = "Black Coucal",
    'WBC' = "White-browed Coucal")

analysis_names <- c(
  'male' = "Male Mo scenario",
  'female' = "Female Mo scenario"
)


# plot the comparative LTRE results
vital_rate_theme <- 
  theme_bw() +
  theme(
    text = element_text(family = "Verdana"),
    legend.position = "none",
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.ticks.length = unit(0.1, "cm"),
    panel.border = element_blank(),
    panel.spacing.x = unit(0.3, "lines"),
    panel.spacing.y = unit(0.7, "lines"),
    strip.background = element_blank()
  )

# theme_set(vital_rate_theme)

LTRE_ASR_plot_background <-
  ggplot2::ggplot(data = LTRE_coucal_ASR,
                  aes(x = parameter, y = contribution, fill = parameter)) +
  annotate("rect", xmin=-Inf, xmax=Inf, ymin=-0.15, ymax=0, alpha=0.6,
           fill= sex_pal2[1]) + 
  annotate("rect", xmin=-Inf, xmax=Inf, ymin=0, ymax=0.15, alpha=0.6,
           fill= sex_pal2[2]) + 
  annotate("text", x = c(2), y = c(-0.14),
           label = c("\u2640"), size = 4, family = "Menlo",
           vjust = c(0), hjust = c(0.5), colour = "grey10") +
  annotate("text", x = c(2), y = c(0.14),
           label = c("\u2642"), size = 4, family = "Menlo",
           vjust = c(1), hjust = c(0.5), colour = "grey10") +
  facet_grid(sex ~ species, 
             labeller = labeller(.cols = species_names, .rows = analysis_names)) +
  vital_rate_theme +
  theme(axis.title.x = element_text(size = 7, vjust = -0.1, colour = "white"),
        axis.text.x  = element_text(size = 7, angle = 45, hjust = 1, colour = "white"),
        axis.ticks.x = element_line(size = 0.2, colour = "white"),
        strip.text.x = element_text(size = 6, colour = "white"),
        axis.title.y = element_text(size = 7, hjust=0.5, vjust = 3.5, colour = "white"),
        axis.text.y  = element_text(size = 6, colour = "white"),
        axis.ticks.y = element_blank(),
        strip.text.y = element_text(size = 6, colour = "white"),
        plot.margin = unit(c(0.2, 0.85,
                             0.78, 0.6), "cm")) +
  
  scale_x_discrete(labels = c("Hatching sex ratio" = expression(italic("\u03C1")),
                              "Nestling survival" = expression(italic(S["n"])),
                              "Groundling survival" = expression(italic(S["g"])),
                              "Fledgling survival" = expression(italic(S["f"])),
                              "Adult survival" = expression(italic(phi["a"])),
                              "Mating system" = expression(italic("h")))) +
  scale_y_continuous(limits = c(-0.15 , 0.15), expand = c(0, 0),
                     breaks = c(-0.15, -0.1, -0.05, 0, 0.05, 0.1, 0.15)) +
  ylab("Absolute contribution to adult sex ratio bias") +
  xlab("Sex bias in parameter")

LTRE_ASR_plot <-
  ggplot2::ggplot() +
  geom_bar(data = LTRE_coucal_ASR,
           aes(x = parameter, y = contribution, fill = parameter), 
           color = "black", stat = "identity", size = 0.2) + 
  facet_grid(sex ~ species, 
             labeller = labeller(.cols = species_names, .rows = analysis_names)) +
  vital_rate_theme +
  theme(axis.title.x = element_text(size=7, vjust=-0.1),
        axis.text.x  = element_text(size=7, angle = 45, hjust = 1),
        axis.ticks.x = element_line(size = 0.2, colour = "grey40"),
        strip.text = element_text(size = 6, colour = "black", face = "italic"),
        panel.background = element_rect(fill = "transparent", colour = NA),
        plot.background = element_rect(fill = "transparent", colour = NA),
        axis.title.y = element_text(size=7, hjust=0.5, vjust = 3.5),
        axis.text.y  = element_text(size=6),
        axis.ticks.y = element_line(size = 0.2, colour = "grey40"),
        strip.text.y = element_text(size = 6),
        plot.margin = unit(c(0.2, 0.85, 0.2, 0.6), "cm")) +
  scale_fill_manual(values = cbPalette) +
  scale_y_continuous(limits = c(-0.15 ,0.15), expand = c(0, 0),
                     breaks = c(-0.15, -0.1, -0.05, 0, 0.05, 0.1, 0.15)) +
  scale_x_discrete(labels = c("Hatching sex ratio" = expression(italic("\u03C1")),
                              "Nestling survival" = expression(italic(S["n"])),
                              "Groundling survival" = expression(italic(S["g"])),
                              "Fledgling survival" = expression(italic(S["f"])),
                              "Adult survival" = expression(italic(phi["a"])),
                              "Mating system" = expression(italic("h")))) +
  ylab("Absolute contribution to adult sex ratio bias") +
  xlab("Sex bias in parameter")

# jpeg(filename = "/Users/luketheduke2/ownCloud/coucal_ASR/R_project/products/figures/LTRE_ASR_plot.jpeg",
#      # compression = "none",
#      width = 4.75,
#      height = 3,
#      units = "in",
#      res = 600)

grid::grid.newpage()
pushViewport( viewport( layout = grid.layout( 1 , 1 , widths = unit( 1 , "npc" ) ) ) )
print(LTRE_ASR_plot_background, newpage = FALSE)
print(LTRE_ASR_plot, newpage = FALSE)
grid::popViewport()

# dev.off()
# plot the comparative LTRE results
vital_rate_theme <- 
  theme_bw() +
  theme(
    text = element_text(family = "Verdana"),
    legend.position = "none",
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.ticks.length = unit(0.1, "cm"),
    panel.border = element_blank(),
    panel.spacing.x = unit(0.3, "lines"),
    panel.spacing.y = unit(0.7, "lines"),
    strip.background = element_blank()
  )

LTRE_ASR_plot_background_female <-
  ggplot2::ggplot(data = filter(LTRE_coucal_ASR, sex == "female"),
                  aes(x = parameter, y = contribution, fill = parameter)) +
  annotate("rect", xmin=-Inf, xmax=Inf, ymin=-0.15, ymax=0, alpha=0.6,
           fill= sex_pal2[1]) +
  annotate("rect", xmin=-Inf, xmax=Inf, ymin=0, ymax=0.15, alpha=0.6,
           fill= sex_pal2[2]) +
  annotate("text", x = c(2), y = c(-0.14),
           label = c("\u2640"), size = 4, family = "Menlo",
           vjust = c(0), hjust = c(0.5), colour = "grey10") +
  annotate("text", x = c(2), y = c(0.14),
           label = c("\u2642"), size = 4, family = "Menlo",
           vjust = c(1), hjust = c(0.5), colour = "grey10") +
  facet_grid(. ~ species, 
             labeller = labeller(.cols = species_names)) +
  vital_rate_theme +
  theme(axis.title.x = element_text(size = 7, vjust = -0.1, colour = "white"),
        axis.text.x  = element_text(size = 7, angle = 45, hjust = 1, colour = "white"),
        axis.ticks.x = element_line(size = 0.2, colour = "white"),
        strip.text.x = element_text(size = 6, colour = "white", face = "italic"),
        
        axis.title.y = element_text(size = 7, hjust=0.5, vjust = 3.5, colour = "white"),
        axis.text.y  = element_text(size = 6, colour = "white"),
        axis.ticks.y = element_blank(),
        strip.text.y = element_text(size = 6, colour = "white"),
        plot.margin = unit(c(0.2, 0.85,
                             0.78, 0.6), "cm")) +
  scale_x_discrete(labels = c("Hatching sex ratio" = expression(italic("\u03C1")),
                              "Nestling survival" = expression(italic(S["n"])),
                              "Groundling survival" = expression(italic(S["g"])),
                              "Fledgling survival" = expression(italic(S["f"])),
                              "Adult survival" = expression(italic(phi["a"])),
                              "Mating system" = expression(italic("h")))) +
  scale_y_continuous(limits = c(-0.15 , 0.15), expand = c(0, 0),
                     breaks = c(-0.15, -0.1, -0.05, 0, 0.05, 0.1, 0.15)) +
  ylab("Absolute contribution to\nadult sex ratio bias") +
  xlab("Sex bias in parameter")

LTRE_ASR_plot_female <-
  ggplot2::ggplot() +
  geom_bar(data = filter(LTRE_coucal_ASR, sex == "female"),
           aes(x = parameter, y = contribution, fill = parameter), 
           color = "black", stat = "identity", size = 0.2) + 
  facet_grid(. ~ species, 
             labeller = labeller(.cols = species_names)) +
  vital_rate_theme +
  theme(axis.title.x = element_text(size=7, vjust=-0.1),
        axis.text.x  = element_text(size=7, angle = 45, hjust = 1, colour = "white"),
        axis.ticks.x = element_line(size = 0.2, colour = "grey40"),
        strip.text = element_text(size = 6, colour = "black", face = "italic"),
        panel.background = element_rect(fill = "transparent", colour = NA),
        plot.background = element_rect(fill = "transparent", colour = NA),
        axis.title.y = element_text(size=7, hjust=0.5, vjust = 3.5),
        axis.text.y  = element_text(size=6),
        axis.ticks.y = element_line(size = 0.2, colour = "grey40"),
        strip.text.y = element_text(size = 6),
        plot.margin = unit(c(0.2, 0.85, 0.2, 0.6), "cm")) +
  scale_fill_manual(values = cbPalette) +
  scale_y_continuous(limits = c(-0.15 ,0.15), expand = c(0, 0),
                     breaks = c(-0.15, -0.1, -0.05, 0, 0.05, 0.1, 0.15)) +
  scale_x_discrete(labels = c("Hatching sex ratio" = expression(italic("\u03C1")),
                              "Nestling survival" = expression(italic(S["n"])),
                              "Groundling survival" = expression(italic(S["g"])),
                              "Fledgling survival" = expression(italic(S["f"])),
                              "Adult survival" = expression(italic(phi["a"])),
                              "Mating system" = expression(italic("h")))) +
  ylab("Absolute contribution to\nadult sex ratio bias") +
  xlab("Sex bias in parameter")

# jpeg(filename = "/Users/luketheduke2/ownCloud/coucal_ASR/R_project/products/figures/LTRE_ASR_plot_female.jpeg",
#      # compression = "none",
#      width = 4,
#      height = 1.5,
#      units = "in",
#      res = 600)

grid.newpage()
pushViewport( viewport( layout = grid.layout( 1 , 1 , widths = unit( 1 , "npc" ) ) ) )
print(LTRE_ASR_plot_background_female, newpage = FALSE)
print(LTRE_ASR_plot_female, newpage = FALSE)
grid::popViewport()

# dev.off()
# Determine how much larger the contribution of each vital rates is compared to juvenile survival
# groundling vs nestling:
LTRE_coucal_ASR_prop <- 
  filter(LTRE_coucal_ASR) %>% 
  dplyr::select(-model) %>% 
  mutate(parameter2 = as.character(parameter)) %>% 
  dplyr::rename(parameter1 = parameter)

prop_contribution_table <- 
  LTRE_coucal_ASR_prop %>% 
  tidyr::expand(., parameter1 = parameter1, parameter2 = parameter1) %>%
  left_join(., dplyr::select(LTRE_coucal_ASR_prop, parameter1, contribution, sex, species), 
            by = c("parameter1")) %>% 
  dplyr::rename(contribution1 = contribution) %>% 
  left_join(., dplyr::select(LTRE_coucal_ASR_prop, parameter2, contribution, sex, species), 
            by = c("parameter2", "species", "sex")) %>% 
  dplyr::rename(contribution2 = contribution) %>% 
  mutate(prop_diff = abs(contribution1) / abs(contribution2),
         prop_diff2 = ifelse(contribution1 < 0 & contribution2 < 0, -(contribution1/contribution2),
                         ifelse(contribution1 < 0 & contribution2 > 0, contribution1/contribution2,
                                ifelse(contribution1 > 0 & contribution2 < 0, -(contribution1/contribution2), contribution1/contribution2))),
         parameter1 = factor(parameter1, 
                             levels = c("Hatching sex ratio",
                                        "Nestling survival",
                                        "Groundling survival",
                                        "Fledgling survival",
                                        "Adult survival",
                                        "Mating system")),
         parameter2 = factor(parameter2, 
                             levels = c("Hatching sex ratio",
                                        "Nestling survival",
                                        "Groundling survival",
                                        "Fledgling survival",
                                        "Adult survival",
                                        "Mating system")))

LTRE_heatmap <- 
  prop_contribution_table %>% 
  filter(parameter1 %in% c("Mating system",
                           "Adult survival",
                           "Fledgling survival",
                           "Groundling survival",
                           "Nestling survival",
                           "Hatching sex ratio") & 
                  parameter2 %in% c("Mating system",
                                    "Adult survival",
                                    "Fledgling survival",
                                    "Groundling survival",
                                    "Nestling survival",
                                    "Hatching sex ratio")) %>% 
  ggplot(., 
         aes(parameter1, parameter2, fill = prop_diff2)) + 
  geom_tile() +
  facet_grid(sex ~ species,
             labeller = labeller(.cols = species_names, .rows = analysis_names)) +
  theme(
    axis.text.x  = element_text(size = 7, angle = 45, hjust = 1),
    axis.text.y  = element_text(size = 7),
    axis.title.x  = element_text(size = 7),
    axis.title.y  = element_text(size = 7),
    legend.position = "bottom",
    legend.box = "horizontal",
    legend.text = element_text(size = 7),
    legend.title = element_text(size = 7),
    axis.ticks = element_blank(),
    panel.border = element_rect(color = "grey40", fill = NA),
    
  ) +
  scale_x_discrete(labels = c("Hatching sex ratio" = expression(italic("\u03C1")),
                              "Nestling survival" = expression(italic(S["n"])),
                              "Groundling survival" = expression(italic(S["g"])),
                              "Fledgling survival" = expression(italic(S["f"])),
                              "Adult survival" = expression(italic(phi["a"])),
                              "Mating system" = expression(italic("h")))) +
  scale_y_discrete(labels = c("Hatching sex ratio" = expression(italic("\u03C1")),
                              "Nestling survival" = expression(italic(S["n"])),
                              "Groundling survival" = expression(italic(S["g"])),
                              "Fledgling survival" = expression(italic(S["f"])),
                              "Adult survival" = expression(italic(phi["a"])),
                              "Mating system" = expression(italic("h")))) +
  scale_fill_gradient2(low = sex_pal2[1],
                       mid = "white",
                      high = sex_pal2[2], 
                      name = "Proportional contribution of A compared to B") +#,
  xlab("Parameter A") +
  ylab("Parameter B")

LTRE_heatmap

# ggsave(LTRE_heatmap,
#        filename = "/Users/luketheduke2/ownCloud/coucal_ASR/R_project/products/figures/LTRE_heatmap.jpeg",
#        # compression = "none",
#        width = 4,
#        height = 4,
#        units = "in",
#        dpi = 600)

# LTRE_heatmap_female <- 
#   prop_contribution_table %>% 
#   mutate(prop_diff2 = ifelse(abs(prop_diff2) == 1, 0, prop_diff2)) %>% 
#   filter(sex == "female") %>% 
#   filter(parameter1 %in% c("Mating system",
#                            "Adult survival",
#                            "Fledgling survival",
#                            "Groundling survival",
#                            "Nestling survival",
#                            "Hatching sex ratio") & 
#                   parameter2 %in% c("Mating system",
#                                     "Adult survival",
#                                     "Fledgling survival",
#                                     "Groundling survival",
#                                     "Nestling survival",
#                                     "Hatching sex ratio")) %>% 
#   ggplot(., 
#          aes(parameter1, parameter2, fill = prop_diff2)) + 
#   geom_tile() +
#   facet_grid(. ~ species,
#              labeller = labeller(.cols = species_names)) +
#   theme(
#     axis.text.x  = element_text(size = 7, angle = 45, hjust = 1),
#     axis.text.y  = element_text(size = 7),
#     axis.title.x  = element_text(size = 7),
#     axis.title.y  = element_text(size = 7),
#     legend.position = "bottom",
#     legend.box = "horizontal",
#     legend.text = element_text(size = 5),
#     legend.title = element_text(size = 6),
#     axis.ticks = element_line(size = 0.2, colour = "grey40"),
#     panel.border = element_rect(color = "grey40", fill = NA),
#     plot.margin = unit(c(0.2, 0.85, 0.2,
#                          0.55), "cm"),
#     strip.text = element_blank()
#   ) +
#   scale_x_discrete(labels = c("HSR" = "hatching\nsex ratio",
#                               "Nestling" = "nestling\nsurvival",
#                               "Groundling" = "groundling\nsurvival",
#                               "Fledgling" = "fledgling\nsurvival",
#                               "Adult" = "adult\nsurvival",
#                               "Mating system" = "Mating system")) +
#   scale_y_discrete(labels = c("Hatching sex ratio" = expression(italic("\u03C1")),
#                               "Nestling survival" = expression(italic(S["n"])),
#                               "Groundling survival" = expression(italic(S["g"])),
#                               "Fledgling survival" = expression(italic(S["f"])),
#                               "Adult survival" = expression(italic(phi["a"])),
#                               "Mating system" = expression(italic("h")))) +
#   scale_fill_gradient2(low = sex_pal2[1],
#                        mid = "white",
#                        high = sex_pal2[2],
#                        name = "Proportional contribution of A compared to B") +#,
#   xlab("Sex-bias in parameter A") +
#   ylab("Sex-bias in\nparameter B\n ")

# LTRE_heatmap_female

# ggsave(LTRE_heatmap_female,
#        filename = "/Users/luketheduke2/ownCloud/coucal_ASR/R_project/products/figures/LTRE_heatmap_female.jpeg",
#        # compression = "none",
#        width = 4,
#        height = 2.5,
#        units = "in",
#        dpi = 600)

Supplementary Analysis 1: sex-specific juvenile movement

coucal_wp <-
  read.csv("data/raw/juvies_2013_2020_waypoints.csv",
           stringsAsFactors = FALSE, header = TRUE) %>%

  # clean up a few inconsistencies in the data
  mutate(timestamp = ifelse(ring_ID == "GN 76660" & date_dec == "20140502", "2014-05-02 10:09:00.000", timestamp),
         location.long = ifelse(location.long == 34.68357, 34.08357, location.long)) %>%
  filter(ring_ID != "GN 76726") %>%
  filter(ring_ID != "GN 76728") %>%
  filter(ring_ID != "GN 76732") %>%
  filter(ring_ID != "GN 76733") %>%
  filter(ring_ID != "GN 76777") %>% # both male and female??
  filter(movebank_event.id != "16897537877") %>%
  filter(movebank_event.id != "16897535798") %>%
  filter(movebank_event.id != "16897537758") %>%
  filter(movebank_event.id != "16897538205") %>%
  
  # specify the timezone of the time stamp
  mutate(timestamp = as.POSIXct(timestamp,
                                tz = "Africa/Nairobi"),

         # make the serial number a 4 digit string
         serial_number = as.character(str_pad(serial_number, 4, pad = "0"))) %>%
  
  # make a julian date time variable for plotting
  mutate(julian_date = as.numeric(format(timestamp, "%j"))) %>% 
  
  # remove observations older than 70 days of age
  filter(age_days < 71) %>%

  # consolidate columns of interest
  dplyr::select(species, year, nest.nr, ring_ID, sex, timestamp, julian_date, 
                age_days, location.long, location.lat, species, site, location, 
                comment)

# specify coordinate columns
coordinates(coucal_wp) <- c("location.long","location.lat")

# define spatial projection as WGS84
proj4string(coucal_wp) <- CRS("+init=epsg:4326")
# specify mapview options for viewing
mapviewOptions(basemaps = c("Esri.WorldImagery"),
               layers.control.pos = "topright",
               legend.pos = "bottomleft")

# check out the spatial data
mapview(coucal_wp[coucal_wp$ring_ID %in% "GN 76751",], zcol = "julian_date",
        layer.name = "WBC radio tagging",
        layers.control.pos = "topright")
# make a spatial trajectory object of each animal's encounter history
tr_coucal_wp <- 
  as.ltraj(xy = sp::coordinates(coucal_wp),
           date = coucal_wp$timestamp,
           id = coucal_wp$ring_ID)

# convert the trajectory object back to a dataframe, rename columns, and
# calculate the cumulative distance traveled over the observation
tr_coucal_wp <- 
  ld(tr_coucal_wp) %>% 
  dplyr::rename(distance = dist,
                dispersion = R2n,
                cardinal_angle = abs.angle,
                relative_angle = rel.angle,
                ring_ID = id) %>% 
  group_by(ring_ID) %>% 
  mutate(cum_distance = cumsum(distance))

# extract the temporal summary info for each bird (min age, min date, max age) 
# and merge to the trajectory dataframe
coucal_wp_clean <- 
  as.data.frame(coucal_wp) %>% 
  group_by(species, ring_ID, sex, nest.nr) %>% 
  dplyr::summarise(min_age = min(age_days),
            min_date = min(timestamp),
            max_age = max(age_days)) %>% 
  arrange(desc(max_age)) %>% 
  left_join(., tr_coucal_wp, by = "ring_ID") %>% 
  
  # calcuate the differences in days between the observation and minimum date
  # to get time since first observation, convert from seconds to days, then add
  # to the minimum age value to extract the age at a given observation
  mutate(date_diff = date - min_date,
         year = year(date)) %>% 
  mutate(date_diff = round(((as.numeric(date_diff)/60)/60)/24)) %>% 
  mutate(age = min_age + date_diff) %>% 
  ungroup() %>% 
  
  # remove all NA rows
  filter(!is.na(cum_distance)) %>% 
  
  # calculate the temporal lag between observations
  group_by(ring_ID) %>% 
  mutate(obs_diff = date_diff - lag(date_diff),
         origin_x = x[which.min(date)],
         origin_y = y[which.min(date)])

# for each bird calculate the daily distance from origin
coucal_wp_clean$distance_origin <- 
  sapply(1:nrow(coucal_wp_clean), function(i)
    spDistsN1(as.matrix(coucal_wp_clean[i, c("x", "y")]), 
              as.matrix(coucal_wp_clean[i, c("origin_x", "origin_y")]), 
              longlat = TRUE))

# transform the distance metrics to meters
coucal_wp_clean <- 
  coucal_wp_clean %>% 
  mutate(distance = distance * 100000,
         distance_origin = distance_origin * 1000,
         cum_distance = cum_distance * 100000)

# check the distribution of distances
hist(log(coucal_wp_clean$distance))

range(coucal_wp_clean$distance)
boxplot.stats(coucal_wp_clean$distance)
# extract individuals that mainly had observations spaced 2 days apart
two_day_rings <- 
  coucal_wp_clean %>% 
  dplyr::select(species, ring_ID, sex, date_diff, obs_diff, age) %>% 
  group_by(ring_ID) %>% 
  dplyr::summarise(mean_obs_diff = mean(obs_diff, na.rm = TRUE),
            median_obs_diff = median(obs_diff, na.rm = TRUE)) %>% 
  arrange(desc(median_obs_diff)) %>% 
  filter(median_obs_diff <= 2)

# subset to only include individuals with mainly observations spaced 2 days apart
coucal_wp_clean <- 
  coucal_wp_clean %>% 
  filter(ring_ID %in% two_day_rings$ring_ID) %>% 
  
  # remove questionable observation 18 km 
  filter(distance < 18000)
#### Modeling sex-specific movements ####
# Model 1 ####
# Does daily distance moved increase with age since leaving the nest
# and does this relationship differ between males and females?

# M1 Black Coucals: strong effect of age, no interaction with sex
mod_BC_dist <- 
  lmer(log(distance) ~ poly(age, 2) * sex + (1|ring_ID),
         data = filter(coucal_wp_clean, species == "BC"))

# effect-size plots
coefplot::coefplot(mod_BC_dist)

plot(allEffects(mod_BC_dist))

# extract fitted values
mod_BC_dist_fits <- 
  as.data.frame(effect(term = "poly(age, 2) * sex", mod = mod_BC_dist, 
                       xlevels = list(age = seq(min(coucal_wp_clean[coucal_wp_clean$species == "BC", "age"]),
                                                max(coucal_wp_clean[coucal_wp_clean$species == "BC", "age"]), 1),
                                      sex = c("Female", "Male"))))
mod_BC_dist_fits$species <- "BC"

# M1 White-browed Coucals: strong effect of age, no interaction with sex
mod_WBC_dist <- 
  lmer(log(distance) ~ poly(age, 2) * sex + 
         (1|ring_ID), 
       data = filter(coucal_wp_clean, species == "WBC" & distance > 0))

# effect-size plots
coefplot::coefplot(mod_WBC_dist)

plot(allEffects(mod_WBC_dist))

# extract fitted values
mod_WBC_dist_fits <- 
  as.data.frame(effect(term = "poly(age, 2) * sex", mod = mod_WBC_dist, 
                       xlevels = list(age = seq(min(coucal_wp_clean[coucal_wp_clean$species == "WBC", "age"]),
                                                max(coucal_wp_clean[coucal_wp_clean$species == "WBC", "age"]), 1),
                                      sex = c("Female", "Male"))))
mod_WBC_dist_fits$species <- "WBC"

mod_dist_fits <- 
  bind_rows(mod_BC_dist_fits, mod_WBC_dist_fits)
# first set up plotting parameters
# define the plotting theme to be used in subsequent ggplots
luke_theme <- 
  theme_bw() +
  theme(
    text = element_text(family = "Verdana"),
    legend.title = element_text(size = 16),
    legend.text = element_text(size = 12),
    axis.title.x = element_text(size = 12),
    axis.text.x  = element_text(size = 12), 
    axis.title.y = element_text(size = 16),
    axis.text.y = element_text(size = 12),
    strip.text = element_text(size = 12),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.ticks = element_line(size = 0.5, colour = "grey40"),
    axis.ticks.length = unit(0.2, "cm"),
    panel.border = element_rect(linetype = "solid", colour = "grey"),
    legend.position = c(0.1, 0.9)
  )

# set plotting color palettes
plot_palette_sex <- RColorBrewer::brewer.pal(8, "Dark2")[c(2,1)]
plot_palette_brood <- RColorBrewer::brewer.pal(7, "Blues")[c(3:7)]

# specify the facet labels for each species
species.labs <- c("Black Coucal", "White-browed Coucal")
names(species.labs) <- c("BC", "WBC")
# Model 2 ####
# Does cumulative distance moved increase with age since leaving the nest
# and does this relationship differ between males and females?

# M2 Black Coucals: strong effect of age, no interaction with sex
mod_BC_cum_dist <- 
  lmer(log(cum_distance) ~ poly(age, 2) * sex + (1|ring_ID),
       data = filter(coucal_wp_clean, species == "BC"))

# effect-size plots
coefplot::coefplot(mod_BC_cum_dist)

plot(allEffects(mod_BC_cum_dist))

# extract fitted values
mod_BC_cum_dist_fits <- 
  as.data.frame(effect(term = "poly(age, 2) * sex", mod = mod_BC_cum_dist, 
                       xlevels = list(age = seq(min(coucal_wp_clean[coucal_wp_clean$species == "BC", "age"]),
                                                max(coucal_wp_clean[coucal_wp_clean$species == "BC", "age"]), 1),
                                      sex = c("Female", "Male"))))
mod_BC_cum_dist_fits$species <- "BC"

# M2 White-browed Coucals: strong effect of age, no interaction with sex
mod_WBC_cum_dist <- 
  lmer(log(cum_distance) ~ poly(age, 2) * sex + 
         (1|ring_ID), 
       data = filter(coucal_wp_clean, species == "WBC" & distance > 0))

# effect-size plots
coefplot::coefplot(mod_WBC_cum_dist)

plot(allEffects(mod_WBC_cum_dist))

# extract fitted values
mod_WBC_cum_dist_fits <- 
  as.data.frame(effect(term = "poly(age, 2) * sex", mod = mod_WBC_cum_dist, 
                       xlevels = list(age = seq(min(coucal_wp_clean[coucal_wp_clean$species == "WBC", "age"]),
                                                max(coucal_wp_clean[coucal_wp_clean$species == "WBC", "age"]), 1),
                                      sex = c("Female", "Male"))))
mod_WBC_cum_dist_fits$species <- "WBC"

mod_cum_dist_fits <- 
  bind_rows(mod_BC_cum_dist_fits, mod_WBC_cum_dist_fits)
# Model 3 ####
# Do individuals move further from the nest as they age
# and does this relationship differ between males and females?

# M3 Black Coucals: strong effect of age, no interaction with sex
mod_BC_disp <- 
  lmer(log(distance_origin) ~ poly(age, 2) * sex + (1|ring_ID),
       data = filter(coucal_wp_clean, species == "BC" & distance_origin > 0))

# effect-size plots
coefplot::coefplot(mod_BC_disp)

plot(allEffects(mod_BC_disp))

# extract fitted values
mod_BC_disp_fits <- 
  as.data.frame(effect(term = "poly(age, 2) * sex", mod = mod_BC_disp, 
                       xlevels = list(age = seq(min(coucal_wp_clean[coucal_wp_clean$species == "BC", "age"]),
                                                max(coucal_wp_clean[coucal_wp_clean$species == "BC", "age"]), 1),
                                      sex = c("Female", "Male"))))
mod_BC_disp_fits$species <- "BC"

# M3 White-browed Coucals: strong effect of age, no interaction with sex
mod_WBC_disp <- 
  lmer(log(distance_origin) ~ poly(age, 2) * sex + 
         (1|ring_ID), 
       data = filter(coucal_wp_clean, species == "WBC" & distance_origin > 0))

# effect-size plots
coefplot::coefplot(mod_WBC_disp)

plot(allEffects(mod_WBC_disp))

# extract fitted values
mod_WBC_disp_fits <- 
  as.data.frame(effect(term = "poly(age, 2) * sex", mod = mod_WBC_disp, 
                       xlevels = list(age = seq(min(coucal_wp_clean[coucal_wp_clean$species == "WBC", "age"]),
                                                max(coucal_wp_clean[coucal_wp_clean$species == "WBC", "age"]), 1),
                                      sex = c("Female", "Male"))))
mod_WBC_disp_fits$species <- "WBC"

mod_disp_fits <- 
  bind_rows(mod_BC_disp_fits, mod_WBC_disp_fits)
# M1 plot
ggplot() +
  luke_theme +
  geom_line(data = coucal_wp_clean, 
            aes(x = age, y = log(distance), group = ring_ID, color = sex),
            alpha = 0.2) +
  geom_line(data = mod_dist_fits, aes(x = age, y = fit, color = sex),
            lwd = 0.5) +
  geom_ribbon(data = mod_dist_fits, aes(x = age, ymax = upper, ymin = lower, fill = sex),
              lwd = 1, alpha = 0.25) +
  facet_grid(species ~ ., 
             labeller = labeller(species = species.labs)) +
  ylab(expression(paste("Distance moved between observations (log10)" %+-%  "95% CI", sep = ""))) +
  xlab("Age since leaving nest") +
  scale_color_manual(values = plot_palette_sex,
                     name = "Sex",
                     breaks = c("Female", "Male"),
                     labels = c("Female", "Male")) +
  scale_fill_manual(values = plot_palette_sex,
                     name = "Sex",
                     breaks = c("Female", "Male"),
                     labels = c("Female", "Male"))

# M2 plot
ggplot() +
  luke_theme +
  geom_line(data = coucal_wp_clean, 
            aes(x = age, y = log(cum_distance), group = ring_ID, color = sex),
            alpha = 0.2) +
  geom_line(data = mod_cum_dist_fits, aes(x = age, y = fit, color = sex),
            lwd = 0.5) +
  geom_ribbon(data = mod_cum_dist_fits, aes(x = age, ymax = upper, ymin = lower, fill = sex),
              lwd = 1, alpha = 0.25) +
  facet_grid(species ~ ., 
             labeller = labeller(species = species.labs)) +
  ylab(expression(paste("Cumulative distance moved (log10)" %+-%  "95% CI", sep = ""))) +
  xlab("Age since leaving nest") +
  scale_color_manual(values = plot_palette_sex,
                     name = "Sex",
                     breaks = c("Female", "Male"),
                     labels = c("Female", "Male")) +
  scale_fill_manual(values = plot_palette_sex,
                    name = "Sex",
                    breaks = c("Female", "Male"),
                    labels = c("Female", "Male"))

# M3 plot
ggplot() +
  luke_theme +
  geom_line(data = filter(coucal_wp_clean, distance_origin > 0), 
            aes(x = age, y = log(distance_origin), group = ring_ID, color = sex),
            alpha = 0.2) +
  geom_line(data = mod_disp_fits, aes(x = age, y = fit, color = sex),
            lwd = 0.5) +
  geom_ribbon(data = mod_disp_fits, aes(x = age, ymax = upper, ymin = lower, fill = sex),
              lwd = 1, alpha = 0.25) +
  facet_grid(species ~ ., 
             labeller = labeller(species = species.labs)) +
  ylab(expression(paste("Daily distance from nest (log10)" %+-%  "95% CI", sep = ""))) +
  xlab("Age since leaving nest") +
  scale_color_manual(values = plot_palette_sex,
                     name = "Sex",
                     breaks = c("Female", "Male"),
                     labels = c("Female", "Male")) +
  scale_fill_manual(values = plot_palette_sex,
                    name = "Sex",
                    breaks = c("Female", "Male"),
                    labels = c("Female", "Male"))

Supplementary Analysis 2: weekly sex-specific adult survival

status_dat_ad <- 
  # read raw data
  read.csv("data/raw/Coucal_adults_survival_2001-2019_20200129.csv", 
           header = TRUE, stringsAsFactors = FALSE, na.strings = c("", " ", "NA")) %>%
  
  # rename ring_ID column
  dplyr::rename(ring_ID = Alu,
                ageC2 = days_from_ringing_until_last_observation) %>% 
  
  # specify date strings properly
  mutate(month_ringed = str_sub(date_dec_initially_ringed, start = 5, end = 6),
         day_ringed = str_sub(date_dec_initially_ringed, start = 7, end = 8),
         year_ringed = str_sub(date_dec_initially_ringed, start = 1, end = 4),
         month_last_ob = str_sub(date_dec_last_observed, start = 5, end = 6),
         day_last_ob = str_sub(date_dec_last_observed, start = 7, end = 8),
         year_last_ob = str_sub(date_dec_initially_ringed, start = 1, end = 4)) %>% 
  mutate(date_ringed = as.Date(paste(year_ringed, month_ringed, day_ringed, sep = "-"), 
                               format = "%Y-%m-%d"),
         date_last_ob = as.Date(paste(year_last_ob, month_last_ob, day_last_ob, sep = "-"), 
                                format = "%Y-%m-%d")) %>% 
  
  # select variables of interest
  dplyr::select(species, ring_ID, sex, year, ageC2, ageC, statusC, date_ringed, date_last_ob) %>% 
  
  # remove all white space from data
  mutate(across(everything(), ~str_trim(.x))) %>% 
  mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>% 
  
  # specify empty data as NA
  mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>% 
  
  # exclude all individuals that died in the nest
  # filter(Fledged_status == "yes") %>% 
  
  # classify columns
  mutate(sex = as.factor(sex),
         ageC2 = as.numeric(ageC2)) %>% 
  
  # remove rows with missing sex, age, and status info
  filter(!is.na(sex) & !is.na(ageC2) & !is.na(statusC)) %>% 
  
  # make a unique id for each individual
  mutate(# create the age of entry into the data (all at age 15)
    entry = 0,
    
    # specify the age of death or censoring
    exit = ageC2,
    
    # make the event numeric and specify if 
    # the individual died (1) or was censored (0)
    event = as.numeric(statusC),
    
    sex = ifelse(str_detect(sex, pattern = "[Ff]emale"), "F", 
                 ifelse(str_detect(sex, pattern = "[Mm]ale"), "M", sex))) %>% 
  
  # consolidate to variables of interest
  dplyr::select(species, ring_ID, year, sex, entry, exit, event, date_ringed, date_last_ob)

status_dat_ad_week <- 
  status_dat_ad %>% 
  dplyr::select(species, ring_ID, sex, date_ringed, date_last_ob) %>% 
  pivot_longer(-c(species, ring_ID, sex), names_to = "status", values_to = "date") %>%
  mutate(year_numeric = as.numeric(lubridate::year(date)),
         week_numeric = as.numeric(strftime(date, format = "%V"))) %>% 
  dplyr::select(species, ring_ID, year_numeric, sex, week_numeric)

# import raw csv data into R
BC_detect_dat_A <- 
  read_xlsx("data/raw/All_coucal_waypoints_2001_2019_20200202.xlsx", na = "NA", col_types = "text") %>% 
  dplyr::select(species, ring_ID, sex, year, site, age_status, date_dec) %>% 
  mutate(month = str_sub(date_dec, start = 5, end = 6),
         day = str_sub(date_dec, start = 7, end = 8)) %>% 
  mutate(date = as.Date(paste(year, month, day, sep = "-"), format = "%Y-%m-%d")) %>% 
  mutate(across(c("sex", "site", "age_status"), tolower)) %>%
  mutate(age_status = ifelse(age_status == "adult", "A", ifelse(age_status == "juvenile", "J", "XXX")),
         ring_ID = str_replace_all(string = ring_ID, fixed(" "), "")) %>% 
  mutate(across(c("species", "ring_ID", "sex", "site", "age_status"), as.factor)) %>% 
  mutate(sex = ifelse(str_detect(sex, pattern = "[Ff]emale"), "F",
                      ifelse(str_detect(sex, pattern = "[Mm]ale"), "M", sex)),
         month_year = format(date, "%Y-%m"),
         month_numeric = as.numeric(month),
         year_numeric = as.numeric(year),
         week_numeric = as.numeric(strftime(date, format = "%V"))) %>% 
  filter(species == "BC") %>% 
  filter(age_status == "A") %>% 
  dplyr::select(species, ring_ID, year_numeric, sex, week_numeric) %>%
  bind_rows(., filter(status_dat_ad_week, species == "BC")) %>% 
  filter(!is.na(week_numeric)) %>% 
  distinct() %>% 
  mutate(week_std = round(scale_by(week_numeric ~ year_numeric, ., scale = 0), digits = 0)) %>% 
  mutate(week_std = week_std + abs(min(week_std, na.rm = TRUE))) %>% 
  dplyr::rename(year = year_numeric)

WBC_detect_dat_A <- 
  read_xlsx("data/raw/All_coucal_waypoints_2001_2019_20200202.xlsx", na = "NA", col_types = "text") %>% 
  dplyr::select(species, ring_ID, sex, year, site, age_status, date_dec) %>% 
  mutate(month = str_sub(date_dec, start = 5, end = 6),
         day = str_sub(date_dec, start = 7, end = 8)) %>% 
  mutate(date = as.Date(paste(year, month, day, sep = "-"), format = "%Y-%m-%d")) %>% 
  mutate(across(c("sex", "site", "age_status"), tolower)) %>%
  mutate(age_status = ifelse(age_status == "adult", "A", ifelse(age_status == "juvenile", "J", "XXX")),
         ring_ID = str_replace_all(string = ring_ID, fixed(" "), "")) %>% 
  mutate(across(c("species", "ring_ID", "sex", "site", "age_status"), as.factor)) %>% 
  mutate(sex = ifelse(str_detect(sex, pattern = "[Ff]emale"), "F",
                      ifelse(str_detect(sex, pattern = "[Mm]ale"), "M", sex)),
         month_year = format(date, "%Y-%m"),
         month_numeric = as.numeric(month),
         year_numeric = as.numeric(year),
         week_numeric = as.numeric(strftime(date, format = "%V"))) %>% 
  filter(species == "WBC") %>% 
  filter(age_status == "A") %>% 
  dplyr::select(species, ring_ID, year_numeric, sex, week_numeric) %>%
  bind_rows(., filter(status_dat_ad_week, species == "WBC")) %>% 
  filter(!is.na(week_numeric)) %>% 
  distinct() %>% 
  mutate(week_std = round(scale_by(week_numeric ~ year_numeric, ., scale = 0), digits = 0)) %>% 
  mutate(week_std = week_std + abs(min(week_std, na.rm = TRUE))) %>% 
  dplyr::rename(year = year_numeric)

# use the BaSTA function "CensusToCaptHist()" function to convert long format
# encounter histories of each individual, to wide format with 1's and 0's for 
# each week of encounter for each year
BC_weekly_A_ch <- 
  CensusToCaptHist(ID = BC_detect_dat_A$ring_ID,
                   d = BC_detect_dat_A$week_std,
                   dformat = "%W", timeInt = "%W") %>% 
  dplyr::rename(ring_ID = ID) %>%
  mutate(ID = rownames(.)) %>% 
  left_join(., dplyr::select(BC_detect_dat_A, ring_ID, sex, year), by = "ring_ID") %>% 
  distinct()

WBC_weekly_A_ch <- 
  CensusToCaptHist(ID = WBC_detect_dat_A$ring_ID,
                   d = WBC_detect_dat_A$week_std,
                   dformat = "%W", timeInt = "%W") %>% 
  dplyr::rename(ring_ID = ID) %>%
  mutate(ID = rownames(.)) %>% 
  left_join(., dplyr::select(WBC_detect_dat_A, ring_ID, sex, year), by = "ring_ID") %>% 
  distinct()

# create a two-character string for each encounter and clean the output
BC_weekly_A_ch_Burnham <- 
  sapply(BC_weekly_A_ch[2:26], function(x) paste(x, "0", sep = "")) %>% 
  cbind(BC_weekly_A_ch[c(1,45:length(BC_weekly_A_ch))]) %>% 
  mutate(across(everything(), ~str_replace(string = .x, pattern = "NA0", "00"))) %>% 
  mutate(across(everything(), ~str_trim(.x))) %>% 
  mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>% 
  mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>% 
  distinct() %>% 
  mutate(sex = as.factor(sex))

WBC_weekly_A_ch_Burnham <- 
  sapply(WBC_weekly_A_ch[2:26], function(x) paste(x, "0", sep = "")) %>% 
  cbind(WBC_weekly_A_ch[c(1,54:length(WBC_weekly_A_ch))]) %>% 
  mutate(across(everything(), ~str_replace(string = .x, pattern = "NA0", "00"))) %>% 
  mutate(across(everything(), ~str_trim(.x))) %>% 
  mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>% 
  mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>% 
  distinct() %>% 
  mutate(sex = as.factor(sex))

# Import status data
status_dat_all_ad <- 
  # read raw data
  read.csv("data/raw/Coucal_adults_survival_2001-2019_20200129.csv", 
           header = TRUE, stringsAsFactors = FALSE, na.strings = c("", " ", "NA")) %>%
  
  # rename ring_ID column
  dplyr::rename(ring_ID = Alu,
                ageC2 = days_from_ringing_until_last_observation) %>% 
  
  # specify date strings properly
  mutate(month_ringed = str_sub(date_dec_initially_ringed, start = 5, end = 6),
         day_ringed = str_sub(date_dec_initially_ringed, start = 7, end = 8),
         year_ringed = str_sub(date_dec_initially_ringed, start = 1, end = 4),
         month_last_ob = str_sub(date_dec_last_observed, start = 5, end = 6),
         day_last_ob = str_sub(date_dec_last_observed, start = 7, end = 8),
         year_last_ob = str_sub(date_dec_initially_ringed, start = 1, end = 4)) %>% 
  mutate(date_ringed = as.Date(paste(year_ringed, month_ringed, day_ringed, sep = "-"), 
                               format = "%Y-%m-%d"),
         date_last_ob = as.Date(paste(year_last_ob, month_last_ob, day_last_ob, sep = "-"), 
                                format = "%Y-%m-%d")) %>% 
  
  # select variables of interest
  dplyr::select(species, ring_ID, sex, year, ageC2, ageC, statusC, date_ringed, date_last_ob) %>% 
  
  # remove all white space from data
  mutate(across(everything(), ~str_trim(.x))) %>% 
  mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>% 
  
  # specify empty data as NA
  mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>% 
  
  # classify columns
  mutate(sex = as.factor(sex),
         ageC2 = as.numeric(ageC2)) %>% 
  
  # remove rows with missing sex, age, and status info
  filter(!is.na(sex) & !is.na(ageC2) & !is.na(statusC)) %>% 
  
  # make a unique id for each individual
  mutate(# create the age of entry into the data (all at age 15)
    entry = 0,
    
    # specify the age of death or censoring
    exit = ageC2,
    
    # make the event numeric and specify if 
    # the individual died (1) or was censored (0)
    event = as.numeric(statusC),
    
    sex = ifelse(str_detect(sex, pattern = "[Ff]emale"), "F", 
                 ifelse(str_detect(sex, pattern = "[Mm]ale"), "M", sex))) %>% 
  
  # consolidate to variables of interest
  dplyr::select(species, ring_ID, year, sex, entry, exit, event, date_ringed, date_last_ob) %>% 
  dplyr::select(species, ring_ID, year, sex, event)

# join detection history with status data
BC_detect_status_join <-
  filter(status_dat_all_ad, species == "BC") %>%
  left_join(BC_weekly_A_ch_Burnham, ., by = c("ring_ID", "year", "sex")) %>% 
  dplyr::select(-species) %>% 
  mutate(event = ifelse(is.na(event), 0, event))

WBC_detect_status_join <-
  filter(status_dat_all_ad, species == "WBC") %>%
  left_join(WBC_weekly_A_ch_Burnham, ., by = c("ring_ID", "year", "sex")) %>% 
  dplyr::select(-species) %>% 
  mutate(event = ifelse(is.na(event), 0, event))

# determine the last and first detection for each individual 
max_age_index_BC <- 
  apply(BC_detect_status_join[, c(1:25)], 1, function(x) which(x == "10")) %>% 
  lapply(., function(x) x[which.max(x)]) %>% 
  unlist(.)

min_age_index_BC <- 
  apply(BC_detect_status_join[, c(1:25)], 1, function(x) which(x == "10")) %>% 
  lapply(., function(x) x[which.min(x)]) %>% 
  unlist(.)

max_age_index_WBC <- 
  apply(WBC_detect_status_join[, c(1:25)], 1, function(x) which(x == "10")) %>% 
  lapply(., function(x) x[which.max(x)]) %>% 
  unlist(.)

min_age_index_WBC <- 
  apply(WBC_detect_status_join[, c(1:25)], 1, function(x) which(x == "10")) %>% 
  lapply(., function(x) x[which.min(x)]) %>% 
  unlist(.)

# put "11" at the last detection if the status was a "1" (i.e., dead)
for(i in 1:nrow(BC_detect_status_join)){
  if(BC_detect_status_join$event[i] == "1"){
    BC_detect_status_join[i, max_age_index_BC[i]] <- "11"
  }
}

for(i in 1:nrow(WBC_detect_status_join)){
  if(WBC_detect_status_join$event[i] == "1"){
    WBC_detect_status_join[i, max_age_index_WBC[i]] <- "11"
  }
}

# bind the min and max age info to the capture history data
detect_status_BC <-
  bind_cols(BC_detect_status_join, max_age_index_BC) %>% 
  bind_cols(., min_age_index_BC) %>% 
  dplyr::rename(std_week_last_obs = '...31',
                std_week_first_obs = '...32')

detect_status_WBC <-
  bind_cols(WBC_detect_status_join, max_age_index_WBC) %>% 
  bind_cols(., min_age_index_WBC) %>% 
  dplyr::rename(std_week_last_obs = '...31',
                std_week_first_obs = '...32')

# consolidate capture histories for each sex
Black_Coucal_adult_weekly_Burnham_ch <-
  apply(detect_status_BC[, c(1:25)], 1, paste, collapse = "") %>% 
  as.data.frame() %>% 
  dplyr::rename(ch = '.') %>% 
  mutate(ch = as.character(ch)) %>% 
  dplyr::bind_cols(., detect_status_BC[, c(26:length(detect_status_BC))]) %>% 
  filter(str_detect(ch, "11", negate = TRUE) | str_count(ch, "1") > 2) %>%
  dplyr::select(ch, ring_ID, sex, year) %>%
  drop_na() %>% 
  mutate(year = ifelse(year == "1006", "2006", year)) %>%
  mutate_at(vars(ring_ID, sex, year), as.factor) 

White_browed_Coucal_adult_weekly_Burnham_ch <-
  apply(detect_status_WBC[, c(1:25)], 1, paste, collapse = "") %>% 
  as.data.frame() %>% 
  dplyr::rename(ch = '.') %>% 
  mutate(ch = as.character(ch)) %>% 
  dplyr::bind_cols(., detect_status_WBC[, c(26:length(detect_status_WBC))]) %>% 
  filter(str_detect(ch, "11", negate = TRUE) | str_count(ch, "1") > 2) %>%
  dplyr::select(ch, ring_ID, sex, year) %>%
  drop_na() %>%
  mutate_at(vars(ring_ID, sex, year), as.factor)

BC_adult_ch <-
  Black_Coucal_adult_weekly_Burnham_ch %>%
  mutate(sex = as.factor(sex),
         year = as.factor(year)) %>%
  dplyr::select(ch, sex, year)

WBC_adult_ch <-
  White_browed_Coucal_adult_weekly_Burnham_ch %>%
  mutate(sex = as.factor(sex),
         year = as.factor(year)) %>%
  dplyr::select(ch, sex, year)

BC_adult_ch %>% 
  mutate(status = ifelse(str_detect(ch, "11"), "dead", "censored")) %>% 
  group_by(sex, status) %>% 
  dplyr::summarise(n_ind = n())

WBC_adult_ch %>% 
  mutate(status = ifelse(str_detect(ch, "11"), "dead", "censored")) %>% 
  group_by(sex, status) %>% 
  dplyr::summarise(n_ind = n())

length(levels(BC_adult_ch$year))
length(levels(WBC_adult_ch$year))
# process the adult data as a "Burnham" analysis
BC_coucal_adult.proc <- RMark::process.data(data = BC_adult_ch,
                                                model = "Burnham",
                                                groups = c("sex", "year"))

WBC_coucal_adult.proc <- RMark::process.data(data = WBC_adult_ch,
                                             model = "Burnham",
                                             groups = c("sex", "year"))

# make design matrix
BC_coucal_adult.ddl <- 
  RMark::make.design.data(BC_coucal_adult.proc)

WBC_coucal_adult.ddl <- 
  RMark::make.design.data(WBC_coucal_adult.proc)

burnham_S_analysis_run_pFr_dot <- 
  function(proc_data, design_data){
    
    # sex-specific model
    S.sex = list(formula = ~sex)
    
    # sex- and week-specific model
    S.sex_Time = list(formula = ~sex + Time)
    
    # sex- and week-specific model (quadratic)
    S.sex_Quad = list(formula = ~sex + (I(Time) + I(Time^2)))
    
    # sex- and week-specific model (interaction)
    S.sex_x_Time = list(formula = ~sex * Time)
    
    # sex- and week-specific model (quadratic with interaction)
    S.sex_x_Quad = list(formula = ~sex * (I(Time) + I(Time^2)))
    
    # p (encounter probability):
    # constant model
    p.dot = list(formula = ~1)
    
    # F (site fidelity probability):
    # constant model
    F.dot = list(formula = ~1)
    
    # r (recovery probability)
    # constant model
    r.dot = list(formula = ~1)
    
    # create model list for all a priori models above
    cml <- RMark::create.model.list("Burnham")
    
    # run the models in program MARK
    model.list <-  RMark::mark.wrapper(cml, data = proc_data, 
                                       ddl = design_data, delete = TRUE, 
                                       wrap = FALSE, threads = 1, brief = TRUE,
                                       silent = TRUE, output = FALSE)
    
    # output the model list and sotre the results
    return(model.list)
  }

adult_S_analysis_out_BC <-
  burnham_S_analysis_run_pFr_dot(proc_data = BC_coucal_adult.proc,
                                 design_data = BC_coucal_adult.ddl)

adult_S_analysis_out_WBC <-
  burnham_S_analysis_run_pFr_dot(proc_data = WBC_coucal_adult.proc,
                                 design_data = WBC_coucal_adult.ddl)

BC_sex_x_quad_reals <- 
  adult_S_analysis_out_BC$S.sex_x_Quad.p.dot.r.dot.F.dot$results$real %>% 
  bind_cols(data.frame(str_split_fixed(rownames(.), " ", 
                                       n = 5)), .) %>% 
  mutate(age = as.integer(unlist(str_extract_all(X4,"[0-9]+"))),
         sex = as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) == "F", 
                                "Female", "Male"))) %>% 
  filter(X1 == "S") %>% 
  dplyr::select(sex, age, estimate, se, lcl, ucl) %>% 
  mutate(species = "Black Coucal")

WBC_sex_x_quad_reals <- 
  adult_S_analysis_out_WBC$S.sex_x_Quad.p.dot.r.dot.F.dot$results$real %>% 
  bind_cols(data.frame(str_split_fixed(rownames(.), " ", 
                                       n = 5)), .) %>% 
  mutate(age = as.integer(unlist(str_extract_all(X4,"[0-9]+"))),
         sex = as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) == "F", 
                                "Female", "Male"))) %>% 
  filter(X1 == "S") %>% 
  dplyr::select(sex, age, estimate, se, lcl, ucl) %>% 
  mutate(species = "White-browed Coucal")

BC_sex_reals <- 
  adult_S_analysis_out_BC$S.sex.p.dot.r.dot.F.dot$results$real %>% 
  bind_cols(data.frame(str_split_fixed(rownames(.), " ", 
                                       n = 5)), .) %>% 
  mutate(age = as.integer(unlist(str_extract_all(X4,"[0-9]+"))),
         sex = as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) == "F", 
                                "Female", "Male"))) %>% 
  filter(X1 == "S") %>% 
  dplyr::select(sex, age, estimate, se, lcl, ucl) %>% 
  mutate(species = "Black Coucal")

WBC_sex_reals <- 
  adult_S_analysis_out_WBC$S.sex.p.dot.r.dot.F.dot$results$real %>% 
  bind_cols(data.frame(str_split_fixed(rownames(.), " ", 
                                       n = 5)), .) %>% 
  mutate(age = as.integer(unlist(str_extract_all(X4,"[0-9]+"))),
         sex = as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) == "F", 
                                "Female", "Male"))) %>% 
  filter(X1 == "S") %>% 
  dplyr::select(sex, age, estimate, se, lcl, ucl) %>% 
  mutate(species = "White-browed Coucal")

time_varying_plot <- 
  ggplot(data = bind_rows(BC_sex_x_quad_reals, WBC_sex_x_quad_reals)) +
    geom_line(aes(x = age, y = estimate, group = sex, color = sex), lwd = 0.8) +
    geom_ribbon(aes(x = age, ymin = lcl, ymax = ucl, group = sex, fill = sex), alpha = 0.5) +
    scale_color_manual(values = sex_pal2,
                       labels = c("Female", "Male")) +
    scale_fill_manual(values = sex_pal2,
                      labels = c("Female", "Male")) +
  facet_grid(. ~ species) +
  theme_bw() +
  theme(text = element_text(family = "Verdana"),
        axis.title = element_text(size = 12),
        axis.text  = element_text(size = 10), 
        strip.text = element_text(size = 12),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_line(size = 0.5, colour = "grey40"),
        axis.ticks.length = unit(0.2, "cm"),
        panel.border = element_rect(linetype = "solid", colour = "grey"),
        legend.position = c(0.8, 0.2),
        legend.title = element_blank()) +
  xlab("Week of breeding season (standardized across years)") +
  ylab("Adult weekly survival probability (± 95% CI)")

constant_plot <- 
  ggplot2::ggplot() +
    geom_errorbar(data = bind_rows(BC_sex_reals, WBC_sex_reals),
                  aes(x = sex, ymax = ucl, ymin = lcl, color = sex),
                  alpha = 1, width = 0.05, lwd = 0.5) +
    geom_point(data = bind_rows(BC_sex_reals, WBC_sex_reals),
               aes(x = sex, y = estimate, fill = sex),
               lwd = 3, shape = 21) +
    scale_color_manual(values = sex_pal2,
                       labels = c("Female", "Male")) +
    scale_fill_manual(values = sex_pal2,
                      labels = c("Female", "Male")) +
    facet_grid(. ~ species) +
    theme_bw() +
    theme(
      text = element_text(family = "Verdana"),
      axis.title = element_blank(),
      axis.text  = element_text(size = 10), 
      strip.text = element_text(size = 12),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      axis.ticks = element_line(size = 0.5, colour = "grey40"),
      axis.ticks.length = unit(0.2, "cm"),
      panel.border = element_rect(linetype = "solid", colour = "grey"),
      legend.position = "none") +
    ylab("Adult weekly survival probability (± 95% CI)") +
    scale_y_continuous(limits = c(0, 1))

combo_plot_S <- 
  time_varying_plot + constant_plot + plot_annotation(tag_levels = "A")

# ggsave(combo_plot_S,
#        filename = "products/figures/combo_plot_S_weekly.jpeg",
#        width = 12,
#        height = 6,
#        units = "in",
#        dpi = 300)

Supplementary Analysis 3: Plumage-based sex-specific age structure

age_plumage_dist_data <-
  read_xlsx("data/raw/age_classes_2001-2020_2.xlsx", na = "NA", col_types = "text") %>%
  dplyr::rename_all(~str_replace_all(., "\\s+", "")) %>%
  dplyr::mutate(species = tolower(species)) %>%
  dplyr::filter(str_detect(string = species, pattern = "black")) %>% 
  dplyr::mutate(month = str_sub(`date(yymmdd)`, start = 5, end = 6),
         day = str_sub(`date(yymmdd)`, start = 7, end = 8)) %>% 
  dplyr::mutate(date = as.Date(paste(year, month, day, sep = "-"), format = "%Y-%m-%d")) %>% 
  # remove all white space from data
  dplyr::mutate(dplyr::across(.cols = dplyr::everything(), 
                str_remove_all, pattern = fixed(" "))) %>% 
  dplyr::mutate(species = "BC",
                sex = ifelse(sex == "female", "F", ifelse(sex == "male", "M", "XXX"))) %>% 
  dplyr::mutate(dplyr::across(.cols = dplyr::everything(), 
                              str_replace_all, pattern = fixed("onebarred"), replacement = fixed("barred"))) %>%
  dplyr::mutate(Secondarycoverts = ifelse(SN == "462", "rufous", Secondarycoverts)) %>% 
  dplyr::mutate(Secondarycoverts = ifelse(SN == "409", "barred", Secondarycoverts)) %>% 
  dplyr::mutate(age_LEP = ifelse(str_detect(primarycoverts, "barred") & 
                                   str_detect(primaries, "barred") & 
                                   str_detect(Secondarycoverts, "barred") & 
                                   str_detect(secondaries, "barred"), 1,
                                 
                                 ifelse(str_detect(primaries, "barred") | 
                                        str_detect(secondaries, "barred"), 2,
                                        
                                        ifelse((str_detect(primarycoverts, "barred") |
                                               str_detect(Secondarycoverts, "barred")) &
                                               (str_detect(primaries, "rufous") |
                                               str_detect(secondaries, "rufous")), 2,
                                               
                                               ifelse(str_detect(primarycoverts, "rufous") & 
                                                         str_detect(primaries, "rufous") & 
                                                         str_detect(Secondarycoverts, "rufous") & 
                                                         str_detect(secondaries, "rufous"), 3, 0))))) %>%
  dplyr::mutate(age = ifelse(age == "juvenile", 1, age)) %>%
  dplyr::mutate(CHECK = ifelse(age != age_LEP, "CHECK", "")) %>% 
  dplyr::mutate(age_LEP = age) %>% 
  dplyr::filter(!is.na(age_LEP)) %>% 
  dplyr::filter(year != "2020") %>% 
  dplyr::group_by(year, sex, age_LEP) %>%
  dplyr::summarise(n_ages = n_distinct(Alu)) %>% 
  dplyr::group_by(year, sex) %>% 
  dplyr::mutate(prop = n_ages/sum(n_ages)) %>%
  dplyr::mutate(age_LEP = as.character(age_LEP)) %>% 
  dplyr::mutate(age_LEP = ifelse(age_LEP == "3", "3+", age_LEP)) %>% 
  dplyr::mutate(age_LEP = fct_rev(as.factor(age_LEP)))

age_plumage_dist_data2 <- 
  age_plumage_dist_data %>% 
  dplyr::group_by(sex, age_LEP) %>% 
  dplyr::summarise(mean_n_ages = mean(n_ages)) %>% 
  dplyr::group_by(sex) %>% 
  dplyr::mutate(prop = mean_n_ages/sum(mean_n_ages),
                year = "Grand\naverage") %>% 
  dplyr::bind_rows(age_plumage_dist_data, .) %>% 
  dplyr::mutate(sex = ifelse(sex == "F", "Female", "Male"))

age_structure_plot <- 
  ggplot(age_plumage_dist_data2, 
       aes(fill = sex, alpha = age_LEP, 
           y = prop, x = sex)) + 
  geom_bar(position="stack", stat="identity") +
  facet_wrap(. ~ year) +
  theme_bw() +
  theme(text = element_text(family = "Verdana"),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.ticks.length = unit(0.1, "cm"),
    panel.border = element_blank(),
    panel.spacing = unit(0.3, "lines"),
    strip.background = element_blank(),
    legend.position = "top") +
  scale_fill_manual(values = sex_pal2, name = "Sex", guide = FALSE) +
  scale_alpha_manual(values = c(0.3, 0.65, 1), name = "Age (years)",
                     guide = guide_legend(title.position = "top", title.hjust = 0.5, label.position = "bottom")) +
  ylab("Proportion of breeding adults") +
  xlab("Sex")

# ggsave(age_structure_plot,
#        filename = "products/figures/age_structure_plot.jpeg",
#        width = 6,
#        height = 6,
#        units = "in",
#        dpi = 300)

age_structure_table <- 
  age_plumage_dist_data2 %>% 
  dplyr::filter(str_detect(year, "Grand")) %>% 
  dplyr::ungroup() %>% 
  dplyr::select(sex, age_LEP, mean_n_ages, prop) %>% 
  dplyr::arrange(sex, desc(age_LEP)) %>% 
  gt(rowname_col = "row",
     groupname_col = "sex") %>%
  cols_label(age_LEP = html("<i>Age class</i>"), 
             mean_n_ages = "Annual\nmean no. ind.",
             prop = "Proportion") %>% 
  fmt_number(columns = vars(mean_n_ages, prop),
             decimals = 2,
             use_seps = FALSE) %>% 
  cols_align(align = "right",
             columns = vars(mean_n_ages)) %>%
  cols_align(align = "left",
             columns = vars(prop)) %>%
  tab_options(row_group.font.weight = "bold",
              row_group.background.color = brewer.pal(9,"Greys")[3],
              table.font.size = 12,
              data_row.padding = 3,
              row_group.padding = 4,
              summary_row.padding = 2,
              column_labels.font.size = 14,
              row_group.font.size = 12,
              table.width = pct(60))

# ggsave(age_structure_table,
#        filename = "products/tables/age_structure_table.jpeg",
#        width = 4,
#        height = 4,
#        units = "in",
#        dpi = 300)